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1.0 INTRODUCTION 


This report presents the final results of contract number 

NAS 9-16281 of the NASA Johnson Space Center. The main goal 

of this research contract was the development of 

mathematical models capable of simulating the transient, 

steady state, and faulted performance characteristics of 

various brushless dc machine-PSA (power switching assembly) 

configurations. These systems are intended for possible 
« 

future use as primemovers in EMAs (electromechanical 
actuators) for flight control applications. These 
machine-PSA configurations include wye, delta, and 
open-delta connected systems, see Figures (1.0-1), (1.0-2), 

and (1.0-3), respectively. The research performed under 
this contract was initially broken down into the following 
six tasks: 

1. Development of mathematical models for various 
machine-PSA configurations, 

2. Experimental validation of the mathematical model, 

3. Shorted turn model for failure modes, 

4. Experimental validation of the mathematical model for 
shorted turn-failure modes, 

5. Tradeoff study, and 

6. Documentation of results and methodology. 
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The fourth task was ultimately eliminated due to problems 
in procuring the required equipment. 

The developed actuator model is very general and can 
handle a large variety of machine-PSA configurations. The 
model includes mutual coupling between the armature phase • 
windings as well as the ability to handle nonlinearities in 
both the inductances and capacitances of the system. This 
program automatically generates the state model from the 
network topology without any tedious hand derivations. This 
model has the capability of handling PSAs with or without 
separate choppers. Furthermore the chopping can be either 
hysteresis or PWM controlled as desired. 

This program has already been applied successfully to 
wye, delta, and open delta connected machines as well as 
machines with shorted turns in the armature windings. 

The model was verified against test data obtained from a 
15 HP samarium cobalt permanent magnet brushless dc machine 
system designed and built for use in electric vehicles. The 
correlation between measured and simulated data was 
excellent in all cases; thereby satisfying the second task 
listed above. 

The third task is concerned with the development and 
verification of a shorted turn model. The shorted turn 
model was successfully implemented and applied to the wye 
connected machine mentioned earlier. Due to the problems in 
securing the required test equipment, the parameters in this 
case were assumed to be linear. 
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Because of the lack of adequate test facilities, the 
fourth task was cancelled. 

The tradeoff study was not completed because of the large 
computer resource outlays required for the finite element 
field analysis and associated calculations. The derivation 
of this model and results obtained are given in the 
following chapters. 
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Figure (l.O-l) Wye Configuration 
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Figure (1,0-2) Delta Configuration 


J 





Figure (1.0-3) Open Delta Configuration 


2.0 THE ACTUATOR MODEL 


A generalized hybrid matrix modeling approach, suitable 
for the simulation of electromechanical systems consisting 
of solid state switching devices, capacitors (linear or 
nonlinear), inductors (linear or nonlinear) as well as 
rotating machines with mutual Coupling between all of the 
machine windings, is presented in this chapter. This 
approach was chosen because it facilitates the automatic 
generation of the system state model directly from the 
network topology. 

The systems that were analyzed using this modeling 
approach were divided into two subsystems: the power 

electronics and the rotating machine. A brief description 
of these two subsystems, as well as their impact on the 
overall modeling approach, is given below: 

1. The power electronics was the first subsystem considered 
here. The state-of-the-art functions of power 
electronics subsystems are to control and process the 
flow of power from the dc source to the machine and vice 
versa. The power diodes and transistors are assumed to 
be the basic switching components of the power 
electronics. These diodes and/or transistors take on 
either very small or very large equivalent resistances 
depending upon their ”on"/"off" status. This 
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nonlinearity requires that the state equations be 
updated whenever a diode or transistor changes state 
(status). The large spread of these resistance values 
causes a correspondingly large spread in the time 
constants of the overall system. Therefore, for this 
class of systems, the governing state equations are very 
stiff. 

2. The second subsystem, which requires special attention, 
is the machine. In general, both the machine winding 
inductances and induced emfs are functions of the 
winding currents (state variables) and rotor angle. The 
nonlinearities associated with these parameters are 
typically much less severe than those produced by the 
switching action of the power diodes within the 
rectifier bridge. 

Both of these subsystems are susceptible to failures or 
faults. Therefore additional consideration was given to the 
allowed possibility of faults in the various components, 
such as the diodes, transistors, fuses, and capacitors of the 
power electronics as well as faults within the machine 
itself. Some of these faults may lead to further stiffness 
of the system state equations. 

The intrinsic characteristics of such electromechanical 
systems, during both normal and faulted operation, requires 
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that the modeling approach be designed to handle very stiff, 

nonlinear differential equations. Hence, the . modeling 

approach presented in this report utilizes a commercially 

available implicit integration routine designed specifically 

for stiff differential equations. This is necessary in 

order to accurately simulate the dynamics of the system with 

a minimum expenditure of computer resources. Also the wide 

variety of topologies encountered in these type of machine 

systems, together with the large number of capacitances and 

inductances associated with these systems, makes it 
* 

mandatory that the state model be generated automatically 
from the system topology and component parameters. 
Therefore the modeling approach, presented in this report 
has the following features: 

• 1. The model is generated automatically from the system 

topology and component parameters; 

2. the model is capable of handling nonlinear component 
parameters; and 

3. the solution routines can handle extremely stiff (ill- 
conditioned) differential equations. 

The theoretical background of a modeling approach with 
these three attributes is the subject of this chapter. The 
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modeling approach makes heavy use of network graph theory 
[1,2], and is based upon a hybrid matrix formulation, [2]. 
This approach is presented in the seven sections of this 
chapter. A brief description of each section follows. 

In the first section, the nonlinear components typically 
found in such electromechanical systems are identified and 
described. A suitable network model for each of these 
components is presented. 

In the second section, the basic concepts of network 
graph theory relevant to the chosen modeling approach are 
discussed. The use of network graph theory allows for 
compact notation and facilitates the actual implementation 
of the model on the digital computer. 

In the third section, the hybrid matrix approach for 
automatically generating the state equations is discussed in 
detail. This approach was chosen since it easily handles 
the various nonlinear components as well as mutual coupling 
present in these types of systems. 
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In the fourth section, the output equations (branch 
voltages and currents) are obtained in terms of the hybrid 
matrix representation used to obtain the state model. The 
output equations express all of the branch voltages and 
currents in terms of the state variables and forcing 
functions. 

The fifth section outlines the procedure used to 
determine the operating point "on'V'off " of the diodes at 
each point in time. The accuracy of the overall solution is 
closely related to the degree of stiffness and the accuracy 
in determining the precise switching points of these diodes. 

The procedure for integrating the state variables forward 
in time is. given in the sixth section. The effects of 
stiffness on the solution accuracy and model performance are 
discussed. 

In the last section, the modeling approach is verified by 
comparison of simulated results with test data obtained from 
a 15 hp samarium cobalt based brushless dc machine designed 
and built for vehicle propulsion. 
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2.1 Models for the Nonlinear Components 

Typical EMA systems contain a number of nonlinearities 
which must be dealt with separately. The nonlinearities of 
greatest concern can be catagorized into four groups: 

1. The first group of nonlinearities is associated with the 
switching action of the power diodes, see Figure 
(2.1-1). This switching action causes extreme changes 
in the electrical time constants of the system. This is 
due to the extremely large change in the conduction 
properties of the diodes between their "on" and "off" 
states. The slope of the switch I-V characteristic, 
shown in Figure (2.1-1), is equal to the reciprocal of 
the equivalent diode resistance at a given operating 
point. The . "on" state of the diode is represented by 
the first quadrant in Figure (2.1-1), where the slope is 
large (low equivalent diode resistance). Conversely, 
the "off" state is represented by the third quadrant of 
this figure. The slope of the diode characteristic, in 
this quadrant, is very small and hence can be 
represented by a very large equivalent resistance. 

The second group of nonlinearities is associated with 
the switching action of the power transistors. This 
switching action also causes extreme changes in the time 
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Fig. (2.1-1) Piecewise Linear Diode Model 
Approximating the I-V Characteristics 
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constants of the electromechanical system. Like the 
diodes, this is due to the extremely large change in the 
conduction properties of the transistors between their 
"on" and "off" states. The "on" state of the transistor 
is represented by a low equivalent resistance. 
Conversely, the ’’off” state is represented by a very 
large equivalent resistance. The only difference 
between the diode and transistor models is the way in 
which the equivalent resistance value is determined 
throughout a given simulation run. 

3. The machine parameters (inductances and emfs) of the ac ' 
machine form the second group of nonlinearities within 
the system. These nonlinearities are due to saturation 
effects within the iron portions of the machine. The 
winding inductances and emfs are influenced also by the 
angular displacement between the rotor and stator. 
Details on the numerical calculation of these parameters 
are given in reference (3, 4, 5] 

4. The last group of nonlinearities considered here is due 
to failures which may occur in the various components of 
the system. The faults or failures considered in this 
work consist of the following: 

a) Open fuses due to electrical or mechanical 
reasons; 

b) open capacitors due to mechanical reasons; 
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c) shorted capacitors due to dielectric breakdown; 

d) open circuited diodes due to electrical or 

mechanical reasons; and 

e) short circuited diodes due to electrical or 

mechanical reasons. 

Nonlinear Network Model of the Diodes and Transistors 

The nonlinearity due to the switching action of a diode 
is modeled by a piecewise-linear representation of the 
specific diode's I-V characteristic, as shown in Figure 
(2.1-1). In terms- of an equivalent network component model, 
the diode is modeled as a nonlinear or bi-valued resistance, 
see reference [6, 7], The "on" (or first quadrant) value of 
this resistance is defined as the forward voltage drop of 
the actual diode divided by the forward diode current at 
rated conditions for the given system under study. This 
data is obtained from the actual manufacturer's data sheets 
of the power diodes. A similar model is used for the 
transistor; however, the resistance value (between collector 
and emitter) is assumed not to be polarity sensitive. 
These simplified switch models were chosen because it was 
found from past experience that the excessive computation 
required to include the dynamics during the actual switching 
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process would add little in the way of improved accuracy in 
the predicted system response. This is primarily due to to 
the fact that the percent of the total simulation duration 
which these switching devices spend within the highly 
nonlinear region of their IV-characteristics is small. This 
is true for these types of machine systems primarily because 
of their relatively low frequency application. The modeling 
approach presented in this report can incorporate equivalent 
network models for these switching devices if additional 
accuracy is required. 

The algorithm used to switch from one diode state to the 
other is described in detail in Section 2.5. On the 
otherhand, the control logic used for determining the 
"on ,, / ,, off " state of each transistor is a function of the 
state variables, forcing functions, the rotor angle, the 
machine speed, the mode of operation, and the commanded 
machine torque. 


Machine Model Including Shorted Turns 

The second nonlinear system component modeled was the 
multiwinding machine. The electromagnetic interactions 
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between the various windings was modeled by means of n- 
coupled coils, one for each of the n windings. The model of 
each winding consists of three series connected components; 
the coil resistance, the incremental self and mutual 
inductances, and the open circuit back emf. The inductances 
and emf's are functions of the level of saturation and hence 
this model is nonlinear. Also, these parameters are 
functions of the rotor angle, hence the model is time 
dependent as well. This machine model is general enough to 
be applied to machines similar to the permanent-magnet dc- 
machine systems described in Chapter 1 . For simplicity, 
the multiple armature machine will be treated as having two 
separate armatures. One of these armatures is assumed to be 
normal or healthy (h) while the other is faulted or 
defective (d). The n-coupled coil approach is quite general 
and can easily be extended to any number of armatures and 
fault configurations. 

The spatial relationships among the axes of the seven 
windings of this machine is displayed in Figure. (2.1-2). 
Notice that there are two windings for each of three phases 
(a,b,c) as well as a seperate field winding (f). The coil 
axis or positive MMF direction of each coil is indicated by 
the arrows labeled (a,b,c,f). The angle between the 
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positive MMF of phase A and positive field MMF is denoted by 
an angle theta. The behavior of the field winding is 
assumed to be controlled by the value of the field current, 
which is treated as an external known forcing function and 
hence not a state variable. In the case of a permanent 
magnet machine, this field current would be held constant. 
Based upon these assumptions, the voltage equations for the 
seven winding machine model can be expressed in matrix form, 
as follows: 



( 2 . 1 - 1 ) 


Notice that the effect of the field winding was incorporated 
into the vector of winding flux linkages since the field 
current was assumed to be a known function of time because 
of the permanent magnets. Equation (2.1-1) can be expressed 
in shorthand matrix notation, specifically advantageous for 
partial fault analysis, as follows: 
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( 2 . 1 - 2 ) 




0 




+d/dt| — | 



where, v^ and v d represent the voltage vectors of the 

healthy (h) and faulted’ or defective (d) armature 
phase windings (terminal voltages), respectively; 

and R d represent the resistance matrices of the 
healthy and faulted phase windings, respectively; 

i^ and i. d represent the current vectors of the 

healthy (h) and faulted (d) phase windings, 
respectively; and, 

X^ and X d represent the flux linkage vectors of the 
healthy and faulted phase windings, respectively. 



The flux linkage vectors, X^ and X^, are functions of the 
six armature currents, (i^..., i g ), the field current, i f , 
and the rotor position, 0. Hence, by the chain rule: 


[dX^/dt] = [ 3X h /ai 1 ]*(di 1 /dt] + ... + [ 3X h /3i 6 ] * [di & /dt] 

+ [3X h /3i f ]*[di f /dt] + [ 3X h /30 ] * [ d0/dt ] (2.1-3) 

[dX d /dt] = [ 3X d /3i 1 ] * (di 1 /dt] + ... + [3X d /3i fi ]*[di 6 /dt] 

+ [3X d /3i f ]*[di f /dt] + [ 3X d /30]*[d0/dt] (2.1-4) 


« 
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The first derivative of the rotor angle, 0, with respect 
to time, t, is defined as the radian velocity, u, of the 
machine as follows: 

w = d0/dt (2.1-5) 


The "induced armature emf" cah be obtained from equations 
(2.1-3) and (2.1-4) as follows: 


! ~ h ! 

| t 8 V 8 i f’j 

1 [ ax, /ao ] | 
1 i 

i i 

1 — * [di f /dt] + 

1 1 

| *dj 

U>V al f ] j 

| 

1 I3X./30] | 
1 1 


( 2 . 1 - 6 ) 


The parti als of X h and X d with respect to the currents are, 
by definition, the machine incremental inductances. These 


inductances are given in matrix form as: L^, L hd / -hf' 

L^, L dd / and L^, and are defined as follows: 


Ifch = | S V Si l | 8 V 8i 2 | 8 V 8i 3 | 


£hd - | 8 V 8i 4 | 8 V 8i s | a V 8i s| 


Lhf = jn h /u £ i 


(2.1-7) 


( 2 . 1 - 8 ) 


(2.1-9) 
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— dh = ! 8 V 8i l 


ha = I i>V 8i 4 


&d£ = ! 8 V 81 f 


8 V 8i 2 ! 8i f /8i 3 ! 


V 8i 5 | 8 V 8i sj 


( 2 . 1 - 10 ) 


( 2 . 1 - 11 ) 


( 2 . 1 - 12 ) 


Using equations (2.1-3) through (2.1-12) one can rewrite the 
armature terminal voltage vector, equation (2.1-2), as 
follows: 
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hih j %id j j^h 

I | *d/dt | - - 

-dh j -dd | j— d 


(2.1-13) 


In expanded form, the terminal, voltages of the armature can 
now be written as follows: 


|v,| 
I I 

! R i 

0 

0 

1 0 

0 

o 1 
1 


| e i 

1 v 2 1 
1 1 

10 

1 

R 2 

0 

1 0 

0 

o 1 

I* 2 ! 

| e 2 

1 v 3 1 
1 1 

|0 

j 

0 

R 3 

1 0 

0 

0 I 
1 

I* 3 j 

J e 3 

1 I 

1 





1 

1 1 

* 1 1 + 

1 

|v 4 l 
1 1 

|0 

1 

0 

0 

jv; 

0 

o l 

j^j 

i e 4 

1 1 
M 

|0 

1 

0 

0 

1 0 

R 5 

o 1 

s 1 

\H\ 

I® 5 

1 v 6 1 
1 & l 

|0 

1 

0 

0 

1 0 

0 

R ej 

1 1 
j ±6 j 

j e 6 


’ll 

L 12 

L !3j 

r-l 

►J 

L 15 

Li6 | 

I 1 ! 

’21 

L 22 

L 23 j 

L 24 

L 25 

L 26j 

1*2 

’31 

L 32 

L 33 j 

L 34 

L 35 

L 36 j 

| l3 

*d/dt| — 
] 1 4 

'41 

L 42 

^43 j 

L 44 

L 45 

^46 j 

‘51 

L 52 

L 53 j 

L S4 

L 55 

L 56 j 

! A 5 

’61 

L 62 

L 63 ! 

L 64 

L 65 

L 66 ! 



(2.1-14) 
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Equation (2.1-14) is the governing matrix equation of the 
seven winding machine. This equation corresponds to the 
lumped parameter network model shown in Figure (2.1-2). 

The vectors and matricies employed in Figure (2.1r2) are 
related to this discusion by the following identities: 





\—2 


hh ! -hd 


and 


!-4! 


!-5 


^hd | ^dd | 


|— 6 


(2.1-15) 


(2.1-16) 


(2.1-17) 
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A procedure used to obtain these emf and inductance 
parameters using the finite element method is described in 
detail in references [8, 9, 10]. 


Nonlinear Network Models of the Faulted Components 

One important application of the simulation model 
presented in this work is to analyze the effects of 
component faults in the PSA. These faults may occur in the 
fuses, the snubber capacitors, and in the diode and 

transistor switches. The failures or faults within these 
components may be due to both electrical and/or mechanical 
reasons. 

Nonlinear Component Model of the Fuses 

The nonlinearity of the fuse model results from the 

incorporation of its failure modes mentioned previously. 

The first fuse failure mode considered here results from the 

z 

fuse thermal rating (i T) being exceeded. This is due to 
sustained excessive fuse current. Once a fuse opens and 
clears, the open fuse is modeled as a very large resistance. 
Thus, the fuse model is somewhat similar to the diode model 
in that it also is modeled as a nonlinear resistor. 

The second failure mode, considered here, results from 
mechanical fatigue. It is assumed that this type of failure 
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mode results in an open (large resistance). The only 

difference in these two failure modes is in their 

initiations. The first failure mode is initiated when the 

2 

accumulated energy exceeds the i T rating of the fuse. The 
second type of failure mode is initiated by the USER at a 
specified time within a given simulation run. 

The use of the low or high values of resistance to 
simulate the state of each fuse makes it possible to use 
just one network topology during both the faulted and 
unfaulted modes of operation. This greatly simplifies the 
analysis. 

The unfaulted state of the fuse is modeled as a very low 

2 

resistance. This resistance value as well as the i T rating 
are obtained directly from the fuse data sheets supplied by 
the manufacturer. An example, [10], of such a data sheet, 
in graphical form, can be seen in Figure (2.1-3). 

The data, given in this figure, which is pertinent to the 
fuse model is summarized below: 

1. THE RATED CURRENT specifies the rated current 

that can pass through the fuse indefinitely without 
causing the fuse to open and clear. 





2 . 


2 

The I T rating of the fuse can be obtained from this 
figure by squaring the RMS value of any current, I 2 , 
within the specified ” section for calculating I T" as 
shown in Figure (2.1-3) and multiplying this result by 
its corresponding time, T 2 , as follows: 

I 2 T = (I 2 ) Z * T 2 (2.1-18) 


The Total Simulation Duration (T gD ) is USER specified 

and is defined as the total duration of the computer 

, simulation run. The intersection of T gD with the 

boundary between safe and open operation regions gives 

the value of the threshold current I^g. This current 

represents the minimum value of the fuse current for 

which the fuse energy is accumulated. If this 

2 

accumulated energy exceeds the i T rating then the fuse 
is opened and cleared. Notice for the specified time 
period, T gD , any value of sustained fuse current below 
*THS not result in a fuse open. 


These parameters are incorporated in the fuse network 
component model for a given simulation. The schematic and 
corresponding network component symbols of this model are 
given in Figure (2.1-4). 
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[Joules] 


g Hysteresis Type of Network 
it Fuse Model 




The latching-hysterisis behavior of the fuse model is 

shown also within this figure. The fuse is modeled as a 

piecewise (latching) linear resistor which can take on two 

distinct and very widely spread values. The nominal 

(unfaulted) value of the fuse resistance can 

obtained from the manufacturer’s data sheet. The fuse will 

always have this low resistance value if the energy 

2 

dissipated within the fuse never exceeds its I T rating. 

This safe region of operation is shown in Figure (2.1-4) and 

corresponds to the region bounded by the shaded (FUSE OPEN 

AND CLEARED) region and the two axes (time versus current) 

as illustrated in Figure (2.1-3). Once the fuse exceeds its 
2 

I T rating, the resistance value of its corresponding 
network component model ( R pusg) is wade very high ( R oPEN^' 
to simulate the open circuit condition. 

The second fuse failure mode, that is the mechanical 

opening of a fuse, is USER initiated. This is accomplished 

by setting the fuse resistance, R EUSE / to be equal to R Q p EN 

at some predetermined time during the computer simulation 

run. These two modes for initiating a fuse failure can be 

2 

considered as natural (I T) and forced (mechanical). The 
latching behavior of the fuse model, as shown in Figure 
(2.1-4), is necessary since once the fuse is opened, it will 
remain open for the remainder of the simulation. 
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Nonlinear Component Model of Faulted Capacitors 

The nonlinearity of the capacitor model is strictly due 
to the incorporation of its failure modes. 

These failure modes are: 

1. An open circuit due to mechanical fatigue; and 

2. a short circuit due to dielectric failure. 

The equivalent network model of the fault capacitor 

during the normal, open, and shorted modes of operation, is 
« 

illustrated in Figure (2.1-5). The two resistors (Rq S and 
R c0 ) within this equivalent component model are bi-valued 
piecewise linear resistors. The very large (open) value and 
the very small (shorted) value of these resistances are 
determined by the specific capacitor mode of operation. 

These three modes of operation of the fault capacitors 
(normal, open, or shorted) can be implemented by the various 
combinations of the resistance values for R CQ and R cs as 
illustrated in Figure (2.1-5). 

The failure modes of the fault capacitor are initiated in 
the following manner: 

1. A shorted capacitor (R QC anc * R SC are set sma ^^ 

values) is initiated if the fault capacitor's voltage 
(which is typically a state variable) exceeds its 
breakdown value. This breakdown voltage value, called 
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V^bd' obtainable from the manufacturer's data sheet. 
Even though a breakdown can occur for different voltages 
(ac or dc), it is the responsibility of the engineer 
when incorporating these voltages to use his/her 
judgement in determining a suitable peak breakdown 
voltage, V CBD / to be used in a given simulation. Once 
the breakdown (shorted fault) has occurred, the values 
of R oc and R sc / in Figure (2.1-5) are latched at very 
small resistive values. 

2 .' The open capacitor (R QC and R gc are set to very large 
values) is assumed to occur due to mechanical fatigue. 
This type of fault is initiated by the USER at a 
specified time during the simulation run. Once the open 
failure mode has occurred, the values of R oc and Rg C are 
latched at very large values of resistance for the 
remaining duration of the simulation. 

The initiation of these two failure modes of the 
capacitors can be considered as 1) natural, in the case of 
the short, and 2) forced, in the case of the open. The 
latching behavior of this model is necessary since it is 
assumed that these failures (faults) are permanent during a 
particular simulation run. 
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a) Unfaulted Capacitor 


C 



b) Shorted Capacitor 


C 



c) Open Capacitor 


C 



Figure (2.1-5) Faulted Capacitor Equivalent Network Model 
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2.2 Network Graph Theory Concepts 


The computer-aided analysis of relatively complex 
systems, such as the EMA systems introduced in Chapter 1, 
requires a systematic procedure for formulating the 
governing state equations on the digital computer. This can 
be accomplished by applying generalized graph theory 
concepts, [1,2], and matrix algebra to the lumped-parameter 
network model of the system under study. The economical 
advantage of this approach becomes apparent when it is 
required to model a number of different topologies (or 
systems) using the same simulation program, as is the case 
here with the various types of power conditioners and 
brushless dc-machine systems considered in this report. 

In this section, a brief background on the fundamental 
concepts and definitions of network graph theory are 
presented. For further detail, see references [1,2]. 

A lumped-parameter network model of a system can be 
represented by means of a directed graph. A directed graph 
consists of a set of branches which are interconnected at 
nodes. A branch is a one-port representation of a network 
component. The relationship between the branch voltage and 
branch current is defined by the component model. Typical 
components encountered in EMA systems include resistors, 
capacitors, inductors, diodes, sources, rotating machines. 


- 31 - 



etc. Models of multi-port components are represented by a 
set of branches (one-ports). 


In this work, the positive direction of the current 

through a branch and the positive direction of the voltage 

rise across a branch are defined using the consumer 

notation, as shown in Figure (2.2-1). In the consumer 

notation, positive branch power is defined as power consumed 

within the branch, thus the branch is said to be acting as a 

sink. On the other hand, negative branch power is defined as 
* 

power generated from within the branch, thus the branch is 
said to be acting as a source. 


Branch power is defined as the product of the branch 
voltage and the branch current. For the k-th branch, as 
shown in Figure (2.2-1), the branch power is defined as: 


p b. 


= v. 


* U 


( 2 . 2 - 1 ) 


The branch voltage, v. , is defined as the difference 

D k 

between the node voltages v N and v N , taking the assumed 

m n 

positive direction (from node m to node n in this case) of 


the current, i. , as follows: 
°k 


v, = v M 
b k N m 


- V. 


N. 


n 


The node voltages v N and v N 


( 2 . 2 - 2 ) 

are defined in terms of a 


m n 

reference (ground) node voltage, v ref . 
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Consider for example, an arbitrary network with a total 
of NB branches and NN nodes. In this network, it is assumed 
that at least two branches are connected to each node. At 
this point, it would be advantageous to define some of the 
basic network graph terms used thoughout the remainder of 
this report. 

Loop A loop is a set of two or more branches which form a 
closed path 

Tree A tree is a subgraph of a network graph which must 
satisfy three conditions. 

1. All NN nodes must be contained in the tree. 

2. All NN nodes must be connected by branches. These 

branches are called twigs. 

3. The twigs may not form any loops 

Cotree The cotree consists of all branches of the original 
network not contained within the chosen tree. The 
cotree branches are called links. The insertion of 
any link into the tree results in a closed loop. 

Cut- Set A cut- set is a "set” of branches which when "cut" 
away from the original network graph will totally 
isolate it into two separate subgraphs. 
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It must be noted that for any network graph there exists 
more than one choice for a tree and its corresponding 
cotree. The particular type of tree which is employed here 
is called a "normal-tree". A normal-tree is defined as a 
tree whose branches are selected in a specified sequence 
according to their component types. 

In particular, the selection of a normal tree requires 
that all the independent voltage sources, E, and capacitors, 
C, within the original network model be chosen with the 
highest priority as twigs. All voltage sources are selected 
first. It is assumed that there are no all voltage source 
loops. Next, all of the capacitors, except for those 
belonging to a C-E loop, are selected as twigs. A C-E loop 
is defined as a loop consisting entirely of capacitors 
and/or independent voltage sources. Only one of the 
capacitors within every C-E loop must be chosen as a link 
(branch of the cotree) and the remaining C-E loop branches 
are chosen as 'twigs (branches of the normal tree). It 
should be noted that only the normal tree capacitor voltages 
are members of the set state variables. A dual situation 
for determining which inductor currents are chosen as state 
variables, is discussed next. 
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The selection of a normal tree requires that no 
independent current sources, J, and as few as possible 
inductors, L, be chosen as twigs. This is to say, the 
selection of the cotree requires that all the independent 
current sources and inductors within the original network be 
chosen with the highest priority as links. All current 
sources are chosen first. It is assumed that there are no 
all-current-source cutsets. Next, all of the inductors, 
except for those belonging to a L-J cutset, are selected as 
links. A L-J cutset is defined as a cutset consisting 
entirely of inductors and/or independent Current sources. 
Only one of the inductors within every L-J cutset must be 
chosen as a twig (branch of the normal tree) and the 
remaining L-J cutset branches are chosen as links (branches 
of the cotree). It should be noted that only the cotree 
inductor currents are members of the set of state variables. 


2.3 Hybrid-Matrix Approach for Automatic Generation 

of the State Model 

Conventional modeling methods based upon frequency-domain 
techniques, when applied to the analysis of electronically 
commutated dc machine systems, often leave much to be 
desired. This is because of the nonlinear components 
associated within these systems, such as the rectifier 
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diodes and the magnetic core of the machine. The 
incorporation of such components into an overall system 
model requires a nonlinear time-domain approach that can 
easily accommodate these component nonlinearities. Also, in 
order to use the many available integration routines, the 
governing differential equations of the actuator network 
model must be explicitly obtained as the normal form state 
equations, [2], as follows: 

p( x ) = f( x,t ) (2.3-1) 

* 

where x is the state vector of the system and f explicitly 
defines the relationship between the first derivative of the 
state vector with respect to time, p( x ), in terms of the 
state vector and the independent variable time, t. It is 
assumed that all the forcing functions are known functions 
of time and/or the state variables, and therefore, are 
incorporated in f of equation (2.3-1). 

In this section, an algorithm which is used to 
automatically generate the normal form state equations from 
the network graph of the system under study, is presented. 
This algorithm is based upon a hybrid matrix approach, 
[2,5], which is well suited for the analysis of nonlinear 
electrical networks. There are basically six steps involved 
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in the formulation of the normal form state model. These 


are : 

Step 1 Translate the electromechanical system into a 
lumped-parameter network model using the component 
models given in Section 2.1. 

Step 2 Choose a "normal tree" for this lumped-parameter 
network model. 

Step 3 From the chosen tree and cotree, obtain a linear 
resistive n-port of this network. Express the 
behavior of this n-port using a hybrid matrix 
formulation. 

Step 4 Determine the operating point of the nonlinear 
resistors. 

Step 5 Solve for the currents of the cotree capacitors and 
the voltages of the tree inductors. 

Step 6 Obtain the normal-form state equations using the 
results of steps three, four, and five. 

These six steps will now be discussed in detail. 

Translate System into Lumped Parameter Network ( Step 1) 

The first step in the derivation of the normal form state 

equations for a system is to choose a suitable lumped 
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parameter network model. In the case of the brushless 
exciter, described in Section 1.2, this consists of: 

1. Replacing the permanent magnet machine with a machine 
model, similar to one given in Figure (2.1-2); 

2. Replacing the components of the power conditioner by 
the components models described in Section (2.1). 

Choose a Normal Tree ( Step 2) 

Once a lumped-parameter network model of the system has 
been obtained, the choice of a normal tree must then be 
considered. A normal tree is a member of a subset of the set 
of all possible trees. The normal tree is chosen using a 
component selection hierarchy. The order or selection 
hierarchy for the twig and link branches must conform to the 
order given in Table (2.3-1). Furthermore, the numbering of 
the branches (twigs and links) must also follow this 
hierarcy. It is extremely important that this hierarchy 
procedure be followed; otherwise, the matrix partitioning 
employed later in this section would be meaningless. 

If the network topology does not change then the branch 
hierarchy and normal tree selection needs to be implemented 
only once per simulation. However, the remaining steps 
relating to component value changes must be implemented each 
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time any change occurs. 


Notice, the branch hierarchy (allowable branch types) 
given in Table (2.3-1) does not include dependent sources. 
-In this report, mutual inductances and capacitances are 
included within the formulation of the state equations. By 
eliminating the dependent sources as allowable components 
within Table (2.3-1), there are no topological restrictions 
placed upon the these mutuals. This is most desirable for 
the types of power-conditioner fed machine systems 
considered in this research. 
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Table (2.3-1) Tree and Cotree Branch Numbering 


ORDER BRANCH TYPE 


VARIABLE 

SUBSCRIPT 


1- st 

2- nd 

3- rd 
' 4-th 

5-th 


6- th 

7- th 

8- th 

9- th 
10-th 


- ( Tree Branches ) 

Independent Voltage Sources 

Capacitors ( linear or nonlinear ) 

Resistors ( nonlinear ) 

1 Inductor ( linear or nonlinear ) 
per ( L-J ) Cutset . 

Resistors ( linear ) 

( Cotree Branches ) 

Resistors ( linear ) 

Independent Current Sources 

Inductors ( linear or nonlinear ) 

Resistors ( nonlinear ) 

1 Capacitor ( linear or nonlinear ) 
per ( C-E ) Loop . 
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Obtain a Linear Resistive n-port from the Chosen Normal 
Tree and Represent the n-port using Hybrid Matrix 
Formulation ( Step 3) 


Once a normal tree has been chosen, the next step is to 
segregate all the linear resistors from the other components 
of the overall network and place them into an n-port sub- 
network as shown in Figure (2.3-1). The motivation for 
forming this n-port is to isolate the linear, time- invariant 
components of the network and incorporate their influence 
upon the other components by using hybrid parameters. Using 
this approach, the states of the nonlinear components become 
boundary conditions for the n-port. Furthermore, the 
relationship between the port voltages and currents are 
defined by means of a linear time invariant hybrid matrix, 
H. This relationship is given by the following matrix 
equation: 



where H is defined in terms of the linear resistors of 
the n-port and is partitioned as follows: 


Itt 

S*. 

H . 

H 

-=t 

— cc 


(2.3-3) 
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TREE BRANDIES 


I 

I 


COTREE BRANCHES 
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The tree, t, and cotree, c, port voltage and current 
vectors of equation (2.3-2) are defined in terms of the port 
vectors shown in Figure (2.3-1) as follows: 



(2.3-4) 


(2.3-5) 



(2.3-6) 
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(2.3-7) 


The H matrix is obtained by applying KCL to the entire 
network in such a way that each twig current is written as a 
linear combination of the link currents as follows: 


< — twigs — > | < — links — > 


= 0 

(2.3-8) 

Notice that the current vectors i R and i G are internal to 
the linear resistive n-port and belong to the tree and 
cotree, respectively. Equation (2.3-8) can be written in a 
more compact form as follows: 


= 0 

(2.3-9) 


< — twigs — > | < — links — > 


D i = 



D i = 


-tt 

-tR 

1 

1 

- 1 - 

StG 

-tc 

-Rt 

-RR 

1 

1 

i 

-RG 

-Rc 


^t 
— R 


±G 


-c 
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Similarly using KVL, one can write the fundamental loop 
equations as follows: 

< — twigs — > | < — links — > 


B v = 


-Gt 

-GR 

1 

1 

-GG 

i 

1 lo 
1 o 
1 o 

5ct 

-cR 

1 

1 

lo 

o 

a 

I 

— cc 


— R 






= 0 


(2.3-10) 


or in compact form as: 

< — twigs — > | < — links — > 


B v = 



= 0 


(2.3-11) 


From network graph theory, [3], one can relate 2 C Q tree and 


—tree as followS: 


B, = -D . 

—tree — cotree 


(2.3-12) 
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Substituting equation (2.3-12) into equation (2.3-10) 
yields: 


B v = 


< — twigs — > < — links — > 


-*£g 


-2r G 



-D^T 0 

I 1 

— tc 

— Rc — cG 

1 

1 

— cc 


-GG 


-Gc 



0 


(2.3-13) 


The relationship between the voltage and current vectors 
of the resistors in the linear resistive n-port network are 
related by Ohms law as follows: 

(2.3-14) 

for the linear tree resistors and 


-R -R -R 


= — G — G 


(2.3-15) 


for the linear cotree resistors. 


Inspection of the hybrid matrix equation of the n-port, 
equation (2.3-2), reveals that the internal variables 
lYG'A( 3 'XR'i£} have been eliminated by incorporating their 
influence, on the port variables, in H. This is accomplished 
by using equations (2.3-8) and (2.3-13) to relate the port 
variables v. and i as shown in equation (2.3-2). 

U C 
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Using the lower partition of equation (2.3-8) and the 
component equations given in equations (2.3-14) and 
(2.3-15 ) , one can write the following: 


— r 1 -R + Erg g 2g * S Rc i c =0 


(2.3-16) 


Solving for v^, using the upper partition of equation 
(2.3-13) yields: 


-G “ -tG -t + -RG -R 


Substituting this equation into (2.3-16) gives: 


(2.3-17) 


Zr 1 Vr * 5r G rl I Et G v t * Erg v r ] ♦ d Rc i c 


= 0 

(2.3-18) 


Solving for v^ in terms of the port variables yields: 

Vr = IZr 1 - Erg z' 1 Erg'” 1 '"Erg 2^ e£g - E Rc i c ) 


(2.3-19) 


Now define a matrix Y as follows: 


Y = [Z" 1 + D__ zZ, 1 D^_ ] 
- -R -RG -G — RG J 


(2.3-20) 
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Using this definition of Y, equation (2.3-19) can be 
simplified to the expression: 



[ -5rg & 


— tG -t 



(2.3-21) 


A similar procedure will now be followed to obtain i G in 
terms of the port variables. Using the upper partition of 
equation (2.3-13) and the component equations given in 
equations (2.3-14) and (2.3-15), one can write the 
following: 


--tG -t " -RG -R -R + -G -G - 


(2.3-22) 


Solving for i R using the lower partition of equation (2.3-8) 
yields: 


— R “-RG -G " -Rc -c 


Substituting this equation into (2.3-22) gives: 

t G + Sr G ?r (5 rg i G + e> Rc i c i ♦ ? G ±g = ° 


(2.3-23) 


(2.3-24) 


Solving for i_ in terms of the port variables yields: 


-G [ -G + -RG -R — RG ^ ^-tG -t " -RG -R -Rc -c ] 


(2.3-25) 
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Now define a matrix Z as follows: 


z = 12g + z r e EG i 


(2.3-26) 


Therefore, equation (2.3-25) can be simplified to 


i G = Z^I-D^c v t - ul G Z R D Rc i e 1 


(2.3-27) 


The final step in obtaining H is accomplished by solving 
i. and v in terms of v. and i . To obtain i. in this form, 

— "C — C — u — C — t 

substitute equation (2.3-27) into the top partition of 
equation (2.3-8). This yields the following matrix 
relationship: 



-S tG Z - 1 IBtc - -RG ?R 2r c i c l 


l-2 tG 2' 1 e? G ) v t 
+ !2 tG z- 1 2 r G z r d Rc - D tc ] i c 



. (2.3-28) 

Similarly, to obtain v , substitute equation (2.3-25) into 

c 

the bottom partition of equation (2.3-13). This yields the 
following matrix relationship: 


50 



-c = etc 2t * ®L S' 1 i-e RG So' 1 etc * t -e Rc ±c j 


= ^tc - e| c r 1 5 rg z^ 1 s* g i v t 
+ i-eL s ' 1 e Rc i i c 


(2.3-29) 


Combining equations (2.3-28) and (2.3-29) into one matrix 
equation gives the following relationship: 



tree response port vector 


cotree response port vector 


i-e tG 2' 1 s>t G i 


[ e tG s ' 1 e RG z R i Rc - e tc ] 


[ 2t=- e Rc S'^-Sro So^tc 1 1 


i-eL x ' 1 e Rc i 



tree stimulus port vector 
cotree stimulus port vector 


(2.3-30) 
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Comparison of this equation with (2.3-2) and (2.3-3) reveals 
that this is the desired hybrid matrix equation originally 
defined in compact form in equation (2.3-2). It is now 
obvious that the influences of the internal linear n-port 
resistors manifest themselves as entries within the hybrid 
matrix, H. The entries of the matrix, H , are functions of 
the network topology and linear resistors only. Therefore, 
for a fixed topology, the matrix, H , needs to be calculated 
only once per simulation run. 

At this point, it is necessary to expand the H matrix, 
defined in equation (2.3-30), into a partitioned form which 
relates all of the voltage and current port vectors of the 
linear resistive n-port, as shown in Figure (2.3-1). This is 
given on the next page as equation (2.3-31). 

A number of the sub-matricies of H given in equation 
(2.3-31) are null. Substitution of these null submatricies 
into equation (2.3-31), on the next page, yields equation 
(2.3-32) . 
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! 1 

! j 

= H 1 1 

which expands to 


i “ c ! 

! ic 1 


~ E - 

i -c. ] -R, 

! V ! 

t 

! -l \%s ! ?c ! 

t 

| t | t 

| t | c 

1 c 1 c I C 1 


j^t 

l §E t c t 

j §E t R t 

j §E t L t 

Vc 

j- E t L c 

'r E t G c 

j^t C c 

k* t 

h't 

j- c t R t 

Vt 

v= 

j sc t L = 

|*v. 

h t c c 

jVt 

k«t 

jS R t R t 

j SR t L t 

v= 

j SR t L C 

| §R t G c 

Kt G c 

| §L t E t 


l §L t R t 

hvt 

§ V= 

v= 

v= 

k«. 

| Sj c E t 

jv* 

k R t 

k«* 

jv= 

K e = 

i SJ = G = 

k c c 

| c c 

| i?L c E t 

j §L c C t 

| §L c R t 


j— L J 
| c c 

j§L l 
1 c c 

k g 

1 C C 

kc c 

| c c 

j §G c E t 

i^C C t 

| S °C R t 

| SG C L t 

J— G J 

j c c 

h G = L = 

! Bg c G c 

! Sg c C c 


§C c C t 

| §C c R t 

hc-t 

i^C J 
| c c 

JSc L 
| c c 

k g 

| c c 

k c c 


1 V 
1 E t 

1 

1 c t 

1 Y 

1 R t 

1 — L 
1 L t 

1 ij 
1 J c 

1 At 
1 L c 

] -G 
1 c 

1 ic 

1 C c 


( 2 . 3 - 31 ) 
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( 2 . 3 - 32 ) 
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Additional properties of the H matrix are given below: 


Stt 

= §« 

(2.3-33) 

See 

= -cc 

(2.3-34) 

Stc 

= -Set 

(2.3-35) 


Determine Operating Point of Nonlinear Resistors ( Step 4) 

The hybrid matrix equation (2.3-31) relates the 
constraints on the port voltages (v ,v. ) and currents 

C U 

(i c ,i t ) imposed by the linear resistive portion of the 
overall network, as illustrated in Eigure (2.3-1). The 
topological constraints imposed by choosing a normal tree, 
as defined in section 2.2, gives the hybrid matrix as 
partitioned in equation (2.3-32). This equation does not 
include the effects of the nonlinear components, or of the 
energy storage elements and independent sources. The effect 
of these variables is included by imposing the I-V 
characteristics (constraints) of these elements as boundary 
conditions to the n (two terminal) ports of the linear time- 
invariant resistive n-port sub-network. 


55 



At any given instant of time, it is assumed that the 

state of the network is fixed; that is, the stored energy of 

each energy storage element remains constant. In this work, 

the state varables are the tree capacitor voltages, v~ , 

c t 

and the cotree inductor currents, i-, . The state 

— c 

(energy) of the network is fixed at any given instant of 
time by holding these state variables constant (along with 
the independent sources). In the case of the nonlinear 
capacitors and inductors, it is assumed that their 
incremental component values are functions of the state 
variables. The derivatives of the state variables as well 
as branch variables of the nonlinear resistors are not 
considered ’'frozen" for any given instant of time, in the 
sense that these variables must be determined iteratively 
for that instant. The term "frozen" is used to represent 
variables that are "known" for a given instant of time. 
Specifically, "frozen" variables are the voltages and/or 
currents which: 

1, are simply a known function of time (forcing functions 
or independent sources), 

2. are known as a result of integration (state 
variables), or 



3. are linear combinations of the state variables and 
forcing functions. 


The I-V characteristics of the nonlinear resistors are 

assumed to be independent of the state variables and forcing 

functions. Thus, the values or operating points of these 

nonlinear resistors must simultaneously satisfy their own 
% 

I-V characteristics as well as the constraints imposed by 
the rest of the network. The procedure for obtaining the 
operating points of the nonlinear resistors will now be 
described. 

As previously mentioned, the nonlinear components as well 

as the independent sources are all connected externally to 

the n (two terminal) ports of the linear resistive n-port, 

as shown in Figure (2.3-1). The influence of the rest of the 

network upon the port varables, i„ and v n , of the 

R t G c 

nonlinear resistors can be isolated out of equation (2.3-32) 
as follows: 
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|H R E 
| R t E t 

|Hr c 
| R t C t 

|H R R 
1 R t R t 

|o R l 

| R t L t 
j 

|Vc 

|Vc 

IVc 

j 

i Pr c 1 
1 R t C c 1 

| 

|-G e. 

| 

c. 

1 

!-g r 


1 

I3g j 

| 

j-G L 

j— G G 

1 1 
!2 0C | 

| c t 

| c t 

| c t 

| c t 

| c c 


| c c 

1 c c 1 

■ 

i 

1 w 
1 >1 

1 

1 

1 c t 

! -K 



! - L 

| i G 

t 

1 Ac 1 

1 c 1 

1 -c 

1 t 

i t 

1 c 

1 c 

1 c 


(2.3-36) 


The presence of the four null submatricies within this 

equation indicates that i„ and v- are independent of v r 

R t G c L t 

and i_ . Therefore, equation (2.3-26) can be reduced as 
~ C c 
follows: 


|h r e 

1 R t E t 
j 

IHr c 
I R t c t 
j 

|Vt 

j 

IVc 

IVc 

!Hr g 1 
1 R t G c 1 

|-G e. 

ISg c. 

j— G R 

j 

I^G j 

| 

ISg l 

1 I 

-G G 

| c t 

| C t 

| c t 

I c c 

| c c 

1 c c 1 






t 

1 — E. 

! -v 

! -K 



! ! 

1 t 

i t 

1 t 

1 c 

1 c 

1 c 1 


(2.3-37) 
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At any given instant of time, the state variables, v_ 

-c t 

and i_ , as well as the forcing functions, v„ and 
_L c -E t 

ij , are fixed (frozen). Hence, the operating points of 
c 

the nonlinear resistive port variables can be determined by 
means of an iterative procedure such that both the I-V 
characteristics of these resistors and equation (2.3-37) 
are satisfied. 

The state variables and forcing functions can be 

considered as fixed boundary conditions to the linear 

resistive n-port. This is accomplished by representing 

these frozen variables as independent sources for that fixed 

instant of time. Therefore, the behavior of the network 

shown in Figure (2.3-1), can be represented for that instant 

of time, as shown in Figure (2.3-2). In this figure, the 

frozen independent port variables are contained within the 

outer box while those which can vary are placed outside. 

The independent port variables, v T and i^ , outside 

b t c 

this box are free to vary. The nonlinear resistors are also 
placed outside this box since their values must be 
determined iteratively, for that given instant in time. 

The operating points of the nonlinear resistors are 

independent of v r and i_ as can be seen in equations 
“ L t _L c 

(2.3-36) and (2.3-37). This is due to the fact that the 
state (stored energy) of the tree inductors and the cotree 
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1REE BRANCHES 


I COTREE BRANCHES 



Figure(2 . 3-2 ) Port Variables held Constant (Frozen) 
while Determining the Operating Points of 
the Nonlinear Resistors ,R. and G 

— u — C 







capacitors are linear functions of the fixed (frozen) state 
variables and independent sources, as given by rows four and 
eight of the matrix equation (2.3-32), respectively. These 
equations in reduced form are: 


-L " — L. J + — L. L — L 

t t C C l c c 


(2.3-38) 


- C c =Sc c E t - E t + - C c C t Yc t 


(2.3-39) 


Therefore, the response of the tree inductors and cotree 
capacitors is known and fixed for the instant in time under 
consideration as depicted in Figure (2.3-3) 

Based upon this discussion, equation (2.3-37) can be 
conceptualized as the sum of a constant source vector, due 
to the fixed state variables and forcing functions 
(represented as "frozen” inputs in Figure (2.3-3)), and a 
vector which is a nonlinear function representing the I-V 
characteristics of the nonlinear resistors. This can be 
represented in matrix form as follows: 


(2.3-40) 
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Figure (2.3-3) Frozen or Fixed Port Variables 

While Determining the Operating Points 
of the Nonlinear Resistors, R. and G . 

U C 
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where the source vector is defined as 


1 s 1 

1 t 1 

i^t c t 

!Vc 

j _ 

IHr e 
1 R t E t 

|Hp j | 

1 Vc 1 

! -g. I 

— | 
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and the vector representing the nonlinear behavior of these 

resistors, obtained from equation (2.3-31) and their 
% 

nonlinear I-V characteristics, is defined as follows: 


- R t 

1 

1 

1 

1 

r 

ft 

! % t G 
1 t 1 

^c 

1 
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i S <v= 

C i 



(2.3-42) 


In this work, it is assumed that the I-V characteristics of 
the nonlinear resistors can be explicitly written as 
follows : 


i Y | 
I t | 


r l< *Rt (1)) 

r 2 ( i Rt ( 2 )) 





(2.3-43) 


r n Rt ( i Rt (n Rt )} 
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(2.3-44) 


1 

1 

9n Gc 

( v Gc ^ n Gc^ * j 




where i^(k) is current through the k nonlinear resistor 
with a corresponding branch voltage of v Rt (k) . 


Therefore, in more compact form, the nonlinear equations 

employed in the iterative process to obtain the operating 

points of the nonlinear resistors, (v., ,i_, ) and 

R t R t 

(y G 'Ac )/ 9 iven as follows: 


= 0 


(2.3-45) 

Note that equations (2.3-43) and (2.3-44) represent 
either the voltage or current of each nonlinear resistor 
explicitly in terms of the other variable. In order to 
obtain a matrix equation, similar to (2.3-45) strictly in 
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terms of one variable (i or v) for each nonlinear resistor, 
all of the nonlinear resistors must have I-V characteristics 
that are single valued. The incorporation of diodes within 
the overall system model does not present any problem since 
their I-V characteristics are monotonic and therefore 
possess an inverse. In other words, either variable (i or 
v) can be written explicitly in terms of the other for the 
diode model presented in this work. 

It is assumed that at each fixed instant of time or 
(equivalently for frozen port sources as shown in Figure 
(2.3-3)), that a unique operating point for all of the 
nonlinear resistors can be determined, [ 2 , 11). Such a 
unique operating point, the solution to equation . (2.3-45) , 
is symbolized as follows: 


= h' RG (2.3-47) 


The primes above these variables indicate that they are the 
unique solution of the nonlinear equation (2.3-45). Notice 
that this solution is uniquely determined by the state 
variables and forcing functions at each instant of time. 
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Solve for the Currents of the Cotree Capacitors and the 
Voltages , of the Tree Inductors ( Step 5) 

Once the operating points of the nonlinear resistors, 

(i 0 ,v„ ) and (i_ ,v- ), within the original 
-R t " R t “ G c "^c 

network have been determined as h* RG ; the unknown tree 

inductor voltages, v r , and the unknown cotree capacitor 

_L t 

currents, i_ , can be determined, independent of the 
~c 

nonlinear resistors, for the given instant in time as 

illustrated in Figure (2.3-2). 

* 

It is assumed, in this work, that component equations 
for all inductors and capacitors (including mutuals) can be 
expressed (on an instant by instant basis) in matrix form 
as follows: 


v Ti | 
t 1 

i htt 

*tc 

1 I Ar 

i 1 L t 




1 * d/dt 1 

! 

! -ct 

^cc 

•H| 

c 1 

1 


1 1 c 


(2.3-48) 


1 in 1 
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• Stt 
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1 1 
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1 1 
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•H | 
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o 

c+ 

£cc 

l< 
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1 c 1 
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1 i 

c 



(2.3-49) 
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The entries of these coefficients matrices contain 
incremental values of inductances and capacitances, 
respectively. The arbitrary (m,n) entries for these two 
coefficient matrices are defined as follows: 


3 < vw » > 

u C 

3 < vv> 

c 

3 < > 

3 < in ) 

3 < %< V C ' V C > > 

3( ,J m (v C >> 

t C 


3 < v n > 

> C 

CO 

1 


(2.3-50) 


(2.3-51) 


where X ( i T , i r ) is the total flux linkage of the 
m “ L t -L c 

m-th inductor and q^Vj, ,y c ) is the total charge of 

t c 

the m-th capacitor. The current, i R , represents the n-th 

current within the set of currents , {i T ,i_ }. 

L t L c 

Similarly the voltage v n represents the n-th voltage within 

the set of capacitor voltages , [v^ , v_ } . Finally, 

c c 

the partial derivative, 3( X m )/3( i R ), in equation 
(2.3-50), represents the instantaneous change of the total 
flux linkage, X , with respect to that n-th winding current, 
i n . Similarily, the partial derivative in equation 
(2.3-51) represents the instantaneous rate of change of the 
charge on the m-th capacitor with respect to the voltage 
across the n-th capacitor. 


In the types of electromechanical systems considered 
here, it was assumed that mutual capacitances could be 
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neglected. In that case, the capacitance matrix, in 
equation (2.3-49), is diagonal. 


Using these component equations and the contraints 
imposed on them by the linear resistive n port; the unknown 
link capacitor currents and twig inductor voltages can now 
be determined. The constraints on these port variables are 
obtained from the fourth and eighth rows of equation 
(2.3-32) as follows: 


% 


±L, 


i^t ! 5 c c L c |- c = E t 


, 1 


|° L t c t |Vc IVt ! s Vc 


I Y c I i L I v„ | i, | 

I c t I L c I E t I J c I 

(2.3-52) 

Differentiating both sides of this equation yields the 
following: 


d/dt 




~L. 


| ic c C t |- C c L c |- C c E t |- C c J c 
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IVt |- L t L c |° L t E t | SL t J c | 


d/dt | y c j i T | Vp | i T | 


— L — E. -J 

t | c | t I c 


(2.3-53) 
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Using equations (2.3-48) and (2.3-49) one can write 
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(2.3-54) 

Substitution of equation (2.3-53) into equation (2.3-54) 
yields: 
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(2.3-55) 
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Collecting terms yields: 
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(2.3-56) 


Thus we have expressed i_ and v r in terms of the 

-C c -L t 

derivatives of both the state variables and forcing 
functions. This completes step 5. 


Derivation of the Normal Form State Equations ( Step 6) 

The final step in obtaining the normal form state 

equations uses the results from steps three, four, and five. 

The response port variables, i_, and v T , of the state 

u t ^c 

variables v„ and i T , respectively are defined by 
~ C t _L c 

rows two and six of equation (2.3-32). 
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These two equations can be expressed as follows: 
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(2.3-57) 


The primes above the nonlinear resistor port variables, 

v' and i’ , indicate the that these satisfy 

R t G c 

equation (2.3-45). 

Next, the substitution of equation (2.3-56) into (2.3-57) 
yields the matrix equation on the next page. 
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(2.3-58) 


vectors i_ and v r as follows: 
~ C t _L c 


72 





Substituting for the derivatives of v_ and i r from 

C c L t 

equation (2.3-53) yields: 
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Regrouping equation (2.3-60) and substituting the result for 
the left side of equation (2.3-58) yields: 
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(2.3-61) 


Finally, by performing the necessary matrix algebra on 
equation (2.3-61), the derivative of the state variables can 
be explicitly expressed in terms of the state variables, the 
forcing functions and their derivatives, as well as the 
(iterated upon) nonlinear resistor port variables as shown 
in equation (2.3-62). 

Thus the general form for equation (2.3-62), where x is 
the state variable vector, u, is the forcing function 
vector, and h’ Rr is the previously solved for nonlinear 
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resistors' port vector as given in equation (2.3-47),' can be 
expressed as follows: 


d( x )/dt = A x + Bu + C d( u )/dt + D h' RG 

(2.3-63) 


t Details on how equation (2.3-62) is numerically 
integrated are discussed in Section 2.5. 

Inspection of the expanded form of the state model, 
equation (2.3-62), reveals that the inverse of the matrix, 
is required. In order to numerically evaluate the above 
matricies {A,B,C, and,D} given in equation (2.3-63), the 
matrix must be inverted only once for networks not 

containing any nonlinear capacitances and/or inductances. 
Fortunately, this matrix has some properties that simplify 
the numerical calculations for its inverse for networks 
containing nonlinear inductors and/or capacitors. These 
properties are discussed next. 
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M cl can be obtained from equation (2.3-62) as follows: 
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(2.3-64) 


All of the hybrid sub-matrices within M CL are skew symmetric 
and hence one can use this information to reduce computer 
memory requirements. Incorporating this fact, M CL can be 
expressed as follows: 
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( 2 . 3 - 65 ) 

From matrix algebra, the inverse of M CL should be 
calculated by taking the inverse of M c and individually 
as follows: 



( 2 . 3 - 66 ) 
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This approach greatly reduces both computer memory and 
computation time to form the required inverse of M CL . 

Further simplifications of M CL are possible if the 
capacitors and inductors are bilateral, or 



9 


$ 





(2.3-67) 

Substitution of these identities into equation (2.3-65) 
yields the following: 
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(2.3-68) 


79 



2.4 Automatic Generation of the Output Equations 

At any given instant in time, the values of all branch 
voltages and currents are completely determined by the state 
variables and forcing functions of the system. In matrix 
form, this statement can be expressed as follows: 

• • 

y = E x + F x + P u + £ u + K h ' RG 

(2.4-1) 


where y is the (output) vector of all the branch voltages 
and currents within the overall network. The vectors, x, u, 
and / were previously defined in equation (2.3-63). 

The dot notation represents the first derivative with 


respect to time (i.e. d(x)/dt = x ). The output vector, 
y, can be expanded in terms of the previously defined branch 
variables, equations (2.3-4) through (2.3-7), as shown in 
Table (2.4-1). 


The tree port voltages, y t , are defined in equation 

(2.3-4). The first two vectors of y t are the independent 

voltage sources, v„ , and the twig capacitor voltages, 

E t 

v>. Remember v_, is one of two components of the 
“ C t C t 

state vector, x, and v„ is one of two components of the 

E t 

forcing function vector u. The other two components of u 

and x are the first two elements of i^ , equation (2.3-7), 

respectively. Again, these four components, v_ , i T , 

“ C t _L c 

y E , and ij , are known for any instant in time. 
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Table (2.4-1) Partitioning of Network Branch Variables 
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The next two components of the output equations are 


u uu Iq , which occur as the third entries of 


and 

z - c 


and i , see equations (2.3-4) and (2.3-7), respectively. 

These two vectors can be obtained by using equations 

(2.3-43) through (2.3-45). Since the nonlinear resistors 

are modeled by piecewise linear, continuous, passive and 

strictly monotonically increasing I-V characteristics, the 

inverse of their component equations exist. Therefore 

v e and i^ can be written as follows: 

R t c 



s 

r_1 (i \ > 

(2.4-2) 


- °c 

= 

9 ‘ 1 te '° e > 

(2.4-3) 

*— * 


Since the operating points of the nonlinear resistors are 
fixed at any given instant in time, the nonlinear operators 
r - ^ and g”^ can - be replaced by the constant coefficient 

matrices R _ o „ and G ^ _ , respectively, for that 
R t R t c c 

instant in time. Thus for that instant, equations (2.4-2) 
and (2.4-3) can be combined into one matrix equation which 
define these nonlinear resistors by the linear relationship 
given in the matrix equation (2.4-4) on the next page. 
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(2.4-4) 


Within this matrix equation, R*"i n and G~i _ are 

" R t R t ' G c G c 

both, in general, constant diagonal matrices only for that 
instant in time. 







Substituting equation (2.4-4) into equation (2.3-45) and 

regrouping, one can write the output equations for v D 

R t 

and i G as shown in equation (2.4-5). The corresponding 
c 

output equations for the variables, i D and v~ , are 

R t G c 

obtained by direct substitution of equation (2.4-5) in 
equation (2.4-4). 


The last components of the port stimulus vectors, y t 
and i , are the twig inductor voltages, v r , and the 

c Li t 

link capacitor currents, i„ These variables are defined 
in equation (2.3-56) in terms of the time derivatives of the 
state variables and the forcing functions. 

Since all of the port stimulus variables, y t and 

1 , are known, the remaining response port variables, 

c 

i. and v , can be obtained directly from equation 

u C 

(2.3-32) . 


83 - 



1 _1 1 1 
- R.R -R G | 

1 

!-c c. 

-1 

I 

1 1 I 

!— C l 1 1 

j 

| t t | t c ! 

1 1 1 

! c t 

1 cell 


1 1-1 I 

-G R - G G 

| c t | c c 1 

|Vt 

1 II 

1 Hr r II 
1 t L c | | 

- 


BR t C t | SR t L c 

1 1 

1 1 

.1*1. 


i IVt 

1 + t 

ll?R J II Z 

I i + 1 

1 

1 

1 

1 


_ 1 

Bc c E t | Bg c L c 

1 1 
1 1 

o 

J 

•H| 

I'l 

1 ISo E 
1 1 G c E t 

i — ii 

| S °c J c 1 1 ±J c 
| c c 1 1 c 

1 

l 

1 

1 



( 2 . 4 - 5 ) 


The only remaining components of y which still need to be 

defined are the branch voltages and currents, v_ , v_ , 

— K — 

i^ , and ig , of the linear resistive n-port. These 
branches are all resistive and internal to the n-port 
subnetwork. 
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The vectors v_. and i_ can be obtained directly in 

terms of the stimulus port vectors, v. and i , using 

equations (2.3-21) and (2.3-25), respectively. The 

associated current and voltage vectors, i R and v^,, are 

defined in terms of v R and i^ by the component equations 

(2.3-14) and (2.3-15), respectively. The component 

coefficient matrices of these linear resistors are diagonal 

and hence their inverses are trivial. This completes the 

process for obtaining the output equations of all the branch 

voltages and currents, y, as defined in equations (2.4-1) 
* 

and (2.4-2). 


2.5 Determination of Network Operating Point 

The state model for many electromechanical systems 
derived in Section (2.3) is nonlinear due to the presence of 
the nonlinear machine parameters (inductances) and solid 
state switching devices (represented by nonlinear 
resistances). This model also has the capability of 
handling nonlinear capacitors; however, this nonlinearity is 
often neglected in the ’ simulation of typical 
electromechanical machine systems. 

The nonlinear inductances are assumed to be functions of 

the state variables, i T , as defined in equation (2.3-50). 

~ L c 

Since the state variables are fixed or frozen at any instant 
of time, the inductances, therefore, are known and fixed at 
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that particular instant of time. The same holds true for 
the nonlinear capacitances. 

The nonlinear resistances, however, must be determined 
iteratively, at each given point in time such that both the 
network constraints and the I-V characteristics of these 
resistors are simultaneously satisfied. 

The nonlinear resistances considered in this work are 
piecewise linear as described earlier in Section (2.1). In 
order to realistically model the switching action of the 
solid state components, a large variation in the slope of 
the I-V characteristics between the "on" and "off" states of 
these switches is required. This slope ratio is of the 
order of six or more magnitudes. Consequently, the 
associated system time constants of the natural response 
also vary by approximately the same orders of magnitude. 
This results in a state model in which the eigenvalues of 
the matrix. A, equation (2.3-63), are very widely separated. 
Such systems are inherently very difficult to solve 
numerically on the digital computer. This problem is 
referred to as the "time constant or stiffness" problem. 

The larger the "off" to "on" ratio of the resistance 
values, the stiffer the governing set of differential 
equations (state model) becomes. Therefore, it is very 
important that the chosen integration routine be able to 
handle this class of differential equations. 
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Furthermore, it will be showh that stiffness affects not 
only the integration process but also the calculation of the 
diode operating points. This is due to the fact that errors 
in the calculation of the diode switching times results in 
entremely large voltage and current transients due to the 
inductive and capacitive nature of the systems under study. 
It will be shown that the "stiffer" the system, the larger 
the magnitude of these fictitious transients. 

An algorithm for determining the operating points of 
these nonlinear resistances will now be discussed. The 
effects of stiffness on the performance of the algorithm 
will also be examined. This algorithm can be broken down 
into the following steps: 

Step 1 Given an updated set of state variables and forcing 
functions at time, t=Tl, determine the present 
values of the nonlinear inductances and 
capacitances. 

Step 2 Solve for the operating points of the nonlinear 
resistances, at t=Tl, by iterating upon equation 
(2.3-45). 

Step 3 Update the matrices given in equation (2.3-62) and 
proceed to integrate the updated state model forward 
in time to t=T2, using an appropriate integration 
routine . 
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Since the solid state switching devices are modeled as 
piecewise linear resistances. Step 2 requires special 
attention at the zero crossing point of the I-V 
characteristic, as illustrated in Figure (2.1-1). 

The most critical step in this algorithm involves the 
precise determination of the all diode switching points 
(times). In previous work, [12-14], these switching points 
were not precisely determined but were merely constrained to 
occur at the chosen solution points. This earlier approach 
was adequate at the time because of the relatively small 
integration steps sizes (5 vs) and because of the relatively 
mild degree of stiffness (R QFF /R 0N > 10000) used by these 
investigators. It was found here, however, that this 
approach for changing the diode states presents serious 
problems for increasing values of "off" resistance and 
increasing integration step sizes. Also, it was found that 
the previous approach, [6-8], required excessive computer 
execution time due to the constraint of smaller integration 
step sizes, required by the explicit integration routine 
employed, for increasing degrees of stiffness. In fact, 
ratios of R 0FF /R QN greater than 10000 resulted in numerical 
problems when using this previous approach. 
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The impact of constraining 'the switching times of the 

diodes to occur precisely at the preceeding integration 

point, which may be different than the true switching times, 

results in the appearance of very large voltage transients 

or spikes within the given numerical simulation run. This 

phenomena can best be understood by examination of the 

example given in Figure (2.5-1). The example consists of a 

series connected voltage source, diode, and inductor in 

which the diode is undergoing a change in state from "on" to 

"off". It is assumed that this change in state occurs 

between the integration points t^ and t^ +1 at a time t=t s . 

Assuming an initial inductor current of 1 ampere at t=t k and 

a sinusoidal source voltage, v_, the diode current is 

s 

assumed to reach a value of -1 ampere at t=t^ + ^ . If the 

diode is constrained to switch at t=t k+1 rather than at the 

correct time, t=t , then a false voltage transient of almost 

s 

10KV is generated across the inductor assuming that v is 

s 

negligible in comparison. In fact, inspection of the 
equation for the inductor voltage, 

v L = L(di/dt) = v s - iRp (2.5-1) 

reveals that the magnitude of this voltage transient is 
approximately a linear function of the "off" resistance. 
Since the magnitude of the series inductor's current is 
constant at any given instant, the change of the equivalent 
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diode resistance at time t=t£ + ^ results in a large increase 
of the diode voltage. This large voltage reflects the 
drastic drop in inductor current as illustrated in Figure 
(2.5-1). The relationship of the integration time step, 
t k+1 -t k , to the voltage spike is not as defined at the "off" 
resistance value of the diode. But, the larger the time 
step, the more negative the inductive current can become, in 
this example, befor the switch in status occurs. This 
larger (more negative) inductor current is directly 
proportional to the voltage spike, as given in equation 
(2.5-1). A dual argument can be made for currents spikes if 
capacitors are present in the system, particlarity if they 
are in parallel with the diode. 

To eliminate these artificial voltage and current 
transients, an improved diode switching algorithm was 
developed. This routine determines all of the diode 
switching times to within a USER specified tolerance. The 
diode status "on"/"off" is only changed at these points in 
time, thereby eliminating the artifical voltage and current 
transients which would have occured using the previous 
approaches. Details of this diode switching algorithm will 


be discussed next. 



+ 




DIODE OFF STATE 
v D < ° 

Bp = lOKft 

V L " V S + 

- v s + 10KV 

- 10KV 

for t “ t X + l 

DIODE ON STATE 
v D> ° 

Ejj = O.Ol 

V L " V S ‘ i D R D 

“ v s + 0.01 
" v s 

for tm i+l 


Figure (2.5-1) Voltage Transient Due to Time Error in Diode 

Switching 


91 - 



Diode Switching Algorithm 


The diode switching algorithm was designed to calculate 
the zero crossing time for each diode. In general these 
switching times may occur in between the discrete 
integration points. Furthermore, this algorithm also 
updates the state model everytime such a change in status is 
detected. The model updating is only done at the calculated 
time of the zero crossing. The basic assumptions upon which 
this algorithm is based are listed below: 

1. The diode components are modeled as piecewise linear 
resistors; 

2. the time interval, At, between integration points is 
small enough so that fast transients are not 
overstepped; 

3. the change of status occurs only at the zero crossings 
of the diodes (within a user specified tolerance); and 

4. a given diode will remain in the "on" state as long as 
its branch voltage is positive and will remain in the 
"off" state as long as its branch voltage is negative. 
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The second assumption is stated to make the reader aware 
of the possibility that this algorithm could miss- a rapid 
change of sign within a diode's voltage or current since a 
polarity comparison is done only at the beginning of the 
time interval. At, and at the end; t fc and t^ + ^ , 
respectively. This blindness is illustrated in Figure 
(2.5-2). Note the polarity check shows that the product of 
the diode current values, at times t^ and is positive, 
thus falsely indicating there was no zero crossing. Since a 
change in polarity would indicate a zero crossing for the 
diode, this particular example illustrates the possibility 
of not detecting such a zero crossing if At is too large. 
Thus, assumption 2 reveals the need to reduce the time 
interval. At, if one desires to track rapid switching 
transients accurately. 


This type of algorithm can be extended to other switching 
devices, such as a silicon controlled rectifiers (SCRs) and 
triacs (back to back SCRs) with some minor modifications. 
These two devices also require the accurate determination of 
their zero crossing times in the same manner as the 
conventional diode, in order to prevent erroneous voltage 
spikes. Details of the necessary logic required for 
modeling such devices will not be discussed here. Also, 
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Figure (2.5-2) Blindness of the Switching Algorithm 
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this algorithm is not limited to one switch point simply at 
the origin. 


The essential features of this algorithm will now be 
discussed. If a change in diode current polarity is 
detected between times t^ and t^.,.-^ for any diode, then 
there exists at least one diode zero crossing within At. 
Once this zero crossing has been detected, the zero crossing 
algorithm (diode switching algorithm) is called. This 
routine is not exited until all of the diodes which cross 
zero within this time interval are identified and dealt with 
properly. This routine entails the following: 


1. Estimate each detected diode zero crossing time. 


2. Rank each crossing time from first to last. 

3. Integrate forward in time, using the unswitched 

(previous) diode status at time t^, to the first 

estimated crossing time, t . 

s 
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4. Check the accuracy of the estimated zero crossing time 
by using the calculated diode currents to determine 
the diode voltage at that point in time. 

5. Change status of the diodes experiencing the zero 

crossing at this time,t , if the diode voltage (and 

s 

current) spikes are within a USER specified tolerance, 
then update the state model. 

6. If t occurs prior to the actual zero crossing time 
consider this estimate to be t^, Figure (2.5-2), and 
go to 1. 

7. In general, either integrate forward or move backward 
in order to determine a more accurate crossing time, 
if needed. 

8. Continue this process until the entire integration 
interval. At, is covered. This . could entail 
subdividing At into a number of subintervals to obtain 
all of the zero crossings. 

Notice that the implementation of item 5 requires updating 
the overall state model due to any change in the diode 
resistances associated with any change in diode status. 

The algorithm employed to model the switching action of 
all the diodes is nested inside the main integration loop of 
the overall simulation program. This algorithm continually 
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checks for any change in polarity between successive diode 
voltage or current values. The major integration loop 
supplies all the branch voltages and currents at every At 
time increment througout the entire simulation duration. 
This algorithm does not interrupt the calculation and 
storage of the solution points at their fixed time interval. 
At. 

Using this improved diode switching algorithm, in 

conjunction with an appropriate implicit integration 

routine, it was possible to more accurately solve stiff 

systems where the ratio between the ''off" and "on" values of 

5 8 

the diode resistance ranged between 10 and 10 . In 
addition, there were savings in computer execution time 
since the implicit routine allows much larger time steps. 
Both of these improvements are discussed in the next 
section. 


2.6 Integration and Stiffness Considerations 

The network operating point and the branch voltages and 

currents, for any given instant in time, are dependent only 

upon the state of the energy storage elements and forcing 

functions (v„ ,i r ,v- , and i T ). Hence in order to obtain 
-C t _L c -E t _J c 

a time solution for all of the branch variables, the state 
equations, given in equation (2.3-62), need to be integrated 
forward in time from one instant to the next using an 
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appropriate integration routine. The selection of an 
appropriate integration routine is based upon solution 
accuracy and the corresponding computational costs for 
extremely stiff systems. 

State equations, associated with power switching type 
networks, are extremely stiff as was mentioned earlier. The 
integration of such stiff differential equations poses a 
number of numerical problems. For example, it is well known 
that the general class of explicit integration routines does 
not compare favorably in terms of accuracy and cost, to the 
implicit predictor-corrector routines when applied to these 
types of systems [5,26]. This and other points will be 
illustrated in the first part of this section by means of a 
numerical example. 

The chosen numerical example is the simple stiff network 
presented in Figure (2.6-1). The numerical solutions from 
these routines are compared with the exact (analytical) 
solution of this network. This comparison is conducted for 
various degrees of stiffness in order to determine which 
type of integration routine is appropriate. 

First, the concept of stiffness will be examined in 
detail using the simple network example. Next, the exact 
solution of this network will be compared with those 
obtained from two commercially available integration 
routines. These two routines are available from the 
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Figure (2.6-1) Example Network Illustrating Stiffness 

of State Equations 
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International Mathematical and ' Statistical Library, IMSL, 
[15]. The chosen explicit routine is a fifth and sixth 
order Runge-Kutta-Verner integration algorithm named DVERK, 
[15, 16]. The chosen implicit routine is a variable order 
Adams predictor-corrector method (or Gears Method) named 
DGEAR, [15, 16]. 

The algorithms employed by these two routines are given 
in Table (2.6-1). These algorithms were obtained from 
reference [17]. DVERK is an explicit routine, by 
definition, since the projected state vector for the next 
time instant, y(t k+1 ), is solely determined by information 
calculated from past history. Considering explicit 
integration algorithms as recursive difference equations 
reveals the relationship between the (maximum allowable) 
step size, h=At, and the algorithm's numerical stability [3, 
5, 20]. DGEAR, on the other hand, is not restricted to 

small integration time steps because of its implicit nature. 

_ Specifically, the homogeneous solution to the difference 
equation for an implicit routine (algorithm) is bounded for 
all positive h. Details of these types of routines can be 
obtained in references [2, 17, 18, 19]. The application and 
blackbox performance of the two integration routines under 
study, for the class of initial valued differential 
equations which are considered stiff, is the concern of this 
section. The concept of stiffness and its probable 
manifestation within a network shall be discussed next. 
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Table (2.6-1) Typical Integration Algorithms 


let, dy/dx = f(x,y) normal-form first-order 

differential equations 


and. 


£ i = f <v*i> 
x i = * h 


f(x 1 ,y(x i )) 


Runge-Kutta algorithms have the form : y^ + ^ = y^ + h 0 

(similar to DVERK) 

where, for a fourth order algorithm; 

0 = 1/6 [k^ + 2(l-.707)k 2 +2(l+.707)k 3 + k 4 ] 

» 

and; 

k l = f i' 

k 2 = f(x.+h/2 , y i +hk 1 /2), 

k 3 = f[x ± +.707h , y i +(-l/2+.707)hk 1 +(l-.707)hk 2 l 
k 4 = f [x ± +h , y i -.707hk 2 +(l+.707)hk 3 ] 


Predictor-Corrector algorithm; (similar to DGEAR) 
Predictor: 

y i+1 = y± W24(55f i -59f i _ 1+ 37f i . 2 -9f._ 3 ) 
Corrector: 

y i+l = y i +h /24(9f i+1 - 1 9f i -5f i-1+ f.. 2 ) 
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The example network, along with its defining state 
equations, is shown in Figure (2.6-1). This particular 
example was chosen because it illustrates how the natural 
response of a network can possess two extremely spread time 
constants. The numerical difficulties encountered with the 
integration of such systems is referred to as the "time 
constant" problem. By increasing the resistance value R 2 , 
which is in series with the inductor, L 2 , the step response 
of the example network becomes more and more difficult to 
calculate accurately using numerical integration routines. 
This is due to the fact that the magnitude of the resistance 
value for R 2 is equal to the stiffness value, o of the 
network under study. 

The stiffness value, o, is defined here as the ratio of 
the two resistor values, o=R 2 /R^. Thus, by increasing only 
one of these resistance values, R 2 for example, while 
holding the other constant (R 1 =l ohm), one can directly vary 
the degree of stiffness of the network example under study. 
The exact relationship between the stiffness value, o, 
defined here, and the ratio of the actual time constants of 
the network's natural response shall be discussed next. 

The analytical unit step response for the two loop 
currents, i^(t) and i 2 (t), shown in Figure (2.6-1) is given 
below: 

il (t) = C 1 e Slt + C 2 e Sjt + a + 1/ Sl s 2 (2.6-2) 
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i 2 ( t) 

= C 3 e Slt + C 4 e S2t + 1/ s x s 2 

(2.6-3) 

where 

- 


C 1 = (2s x + 

0 + l)/(s 1 -s 2 ) 

(2.6-4) 

C 2 = (2s x + 

0 + 1)/(S!-S 2 ) 

(2.6-5) 

C 3 = (s x + 

l)/Sl(Si-S 2 ) 

(2.6-6) 

C 4 = (s a + 

1 )/ S 2 ( Sl”S 2 ) 

(2.6-7) 

Si/S 2 = - (o+3 )/2 

± . 5 [o 2 - 2a + 5] * 5 

(2.6-8) 

Furthermore, 

the steady state response due to 

the unit 

step input is 



lit-) 

= (O+lJ/SiSz 

(2.6-9) 

i 2 (“> 

= (l)/s x s 2 

(2.6-10) 


Examination of equation (2.6-8) reveals that for an 
increasing stiffness value, a, the poles or time constants, 
which are symbolized as s x and s 2 , approach -o and -3/2, 
respectively. Also, the root loci for these two poles 
always remain on the negative real axis of the s-plane 
(Laplace transformation) . Notice also that these two poles 
spread apart (system becomes stiffer), as shown in Figure 
(2.6-2), when o is increased from zero to infinity. 

To characterize the concept of stiffness further; one can 
express the stiffness of the network's state equation in 
terms of the maximum ratio of its "time constants" (or 
eigenvalues of the matrix A in Figure (2.6-1)). 
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The "stiffness" of a set of first order differential 
equations can also be associated with the degree of ill- 
conditioning of the A matrix. More precisely, ill- 

conditioning of a matrix is defined in terms of the ratio of 
its largest to smallest eigenvalue magnitudes. The larger 
this ratio, the more ill-conditioned the A matrix. There is 
a condition number which can be assigned to any matrix in 
terms of its largest and smallest eigenvalues that indicates 
the degree of ill-conditioning of the matrix, [7]. This 

number is sometimes defined as follows: 

| | Xlarge | | 

CN = ■ (2.6-11) 

| | Xsmall | | 

where 

||Xlarge|| is the modulus of the largest eigenvalue 
of the A matrix 

and 


| | Xsmall | | is the modulus of the smallest eigenvalue 
of the A matrix 


In the network example presented in this section, there 
are only two eigenvalues and they are real negative numbers. 
Hence, the condition number, as defined in equation 
(2.6-11), for the A matrix given in Figure (2.6-1) can be 
written as follows: 

CN = [sx/sj] = [ (o+3) + (a 2 -2a + 5 ) ’ 5 /( (o+3 )-(o 2 -2o +5 ) * 5 ) 

( 2 . 6 - 12 ) 
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Notice that as a is increased, the condition number, CN, 
of the A matrix also increases. A large condition number 
indicates difficulties in numerical calculations associated 
the A matrix, such as its inverse for example. A good 
estimate for the number of significant digits of accuracy 
when numerically calculating the inverse of the A matrix can 
be related to the condition number approximately, [17], as 
follows: 

NSD = WLENGTH - [ INT(LOG 10 (CN) ) + 1 ] 

where: 

NSD = number of significant digits 

WLENGTH = number of digits (decimal places) stored in a computer 
word used in computation. 

a) Single Precision = 7 digits 

b) Double Precision = 14 digits 

INT (*) = Truncates the fractional part of the real number, 
hence making it an integer. 

CN = the condition number of the matrix, A. 

For example, if the stiffness value, o, is set equal to 
1000, then the CN is equal to 500.75 and the number of 
digits of round-off that can be expected is 3. On the other 
hand, if a is set equal to 10^ then CN is equal to 
0.50000075 x 10 6 and the number of digits of round-off that 
can be expected when taking the inverse of the A matrix is 
6. These two computations illustrate the direct significance 
of the stiffness value, o, with the round-off error expected 
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using the digital computer. Further, in this particular 
example network, INTflog^CN) + 1 represents the number of 
digits round-off to be expected when taking the inverse of 
the A matrix. 

Similarly, the ill conditioning of some matrix M can also 
have pronounced effects when operating numerically on 
systems of the form M x=b; where b is the vector of known 
forcing functions or measured data and x is the response of 
the system (characterized by M) to the data b. These 
pronounced effects, manifested as large changes in the 
response x due to very small variations or errors in the 
data b, may be caused by numerical round-off of the measured 
data for example. This sensitivity to the data is a direct 
result of the wide spread in the magnitudes of the 
eigenvalues of the matrix M. If a perturbation, Ab,of the 
observed data, b, is strictly in the direction of the 
largest eigenvectors of the matrix M, the resulting solution 
vector x + Ax will not change drastically as compared to 
the case if the same perturbation magnitude occurs in the 


components 

of b 

that 

are 

in the 

same 

direction as 

the 

smallest 

eigenvectors 

of 

M. For 

the 

latter 

case. 

the 

resulting 

change 

of 

the 

solution 

vector, x. 

would 

be 


drastically increased, with respect to the first case, by 
approximately the condition number of the M matrix times the 
perturbation of the data, Ab. 
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The condition number, CN, and the stiffness factor, o, 
are closely related. The context in which they are used is 
their only difference. In particular, the condition number 
of the M matrix, CN, reveals valuable information about 
expected round-off when taking the inverse of M. Similarily 
the stiffness factor, o, gives a measure of how an extreme 
spread in the component values affecting the eigenvalues of 
A can effect the numerical integration of the state 
equations (time constant problem). 

' .Thus, special attention must be given to the specific 
integration routine employed when the A matrix within the 
governing state equation becomes ill-conditioned (a large 
ratio of its largest to smallest eigenvalue magnitudes). 
The larger this ratio, the more ill-conditioned the A 
matrix, and consequently the stiffer the governing state 
equations which contain A become. 

The degree of system stiffness has a great impact on the 
efficiency (computer costs per solution) and the accuracy of 
the numerical integration routines. This is because, in 
order to track the fastest time constant in a stable and/or 
accurate sense, small time steps are required throughout the 
entire simulation period. This is the case even though the 
fastest transient may have died out and consequently no 
longer contributes significantly to the system’ s long term 
response associated with the slower time constants. More 
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precisely/ the maximum allowable integration step size is 
determined by the smallest time constant (or eigenvalue of 
A), while the total number of integration steps required to 
reach steady state is strictly determined by the largest 
time constant (or eigenvalue of A). Hence, whenever the 
time constants (eigenvalues of the A matrix) differ by many 
orders of magnitude, many integration routines (specifically 
explicit . ones) suffer from the requirement of excess 
computer time in order to obtain a desired amount of 
simulation time. 

In order to illustrate the effects of stiffness, two 
commercially available IMSL routines were used to obtain 
numerical solutions for the example network of Figure 
(2.6-1). The defining system of first order differential 
equations, for this network, can be set to any desired 
degree of stiffness by varying the value of o. The two IMSL 
routines used in this example are DVERK and DGEAR. DVERK is 
an explicit initial- value differential equation solver which 
employs a 5th and 6th order Runge-Kutta variable step 
method. DGEAR, on the other hand, is an implicit initial 
value differential equation solver using a Predictor- 
Corrector method with variable step size capability. 
DGEAR' s integration step size does not have to be limited to 
guarantee numerical stability. DVERK, on the other hand, is 
limited to a maximum step size in order to guarantee 
numerical stability. This restriction is inherent in all 



explicit integration routines. Even though DVERK has the 

internal capability to automatically maximize its step size 

throughout a given simulation run, it will not exceed an 

internally generated upper limit in order to guarantee 

numerical stability and/or satisfy the programmer's 

specified maximum local error estimate, TOL. Hence, 

internal to DVERK is an error estimator which limits the 

maximum step size in order to satisfy a specified degree of 

accuracy. DGEAR, on the other hand, can exceed DVERK' s 

maximum step size and still satisfy the same accuracy 
% 

requirements. Details of the numerical analysis of both 
methods can be obtained in references [3, 5, 20]. 

The exact solution to the example network is given in 

equations (2.6-2) and (2.6-3). The corresponding numerical 

6 —7 

solutions from DVERK and DGEAR for o=10 , and T0L=10 , are 

presented in Figures (2.6-3) through (2.6-6). The exact 
solution, in each of these cases, is plotted with the 
corresponding numerical solution for comparison. The global 
relative error, which is defined as the accumulated local 
relative error, is presented within these figures as the 
smaller graph. The global relative error for each of the 
figures is defined as follows: 

(Exact i 1 (t) - DGEAR' s i^t)) 

Relative Error (2.6-13) 

Exact i 1 (t) 
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Relative Error =■ 


(2.6-14) 


(Exact i 2 (t) - DGEAR's i 2 (tT)) 

Exact i 2 (t) 

(Exact i 1 (t) - DVERK's i^T)) 

Relative Error = (2.6-15) 

Exact i^(t) 

(Exact i 2 (t) - DVERK's i 2 (T)) 

Relative Error = (2.6-16) 

Exact i^(t) 

. A comparison of the computation times, memory 
* 

requirements, and solution accuracies of the two IMSL 
routines is presented in Table (2.6-2) for three different 
degrees of stiffness. Based upon these results, it can be 
seen that DGEAR, which is designed specifically for stiff 
systems, out performed DVERK in all of the above catagories. 

Specifically, for a stiffness, o=10 6 , and a specified 

_7 

local relative-error tolerance, TOL=10 , it can be seen 
from this table that DVERK could not keep its maximum global 
relative error (MGRE) within this specified tolerance. 
DGEAR, on the other hand, did a much better job of tracking 
the exact solution. In fact, with this value of stiffness, 
the MGRE of DVERK is nearly one thousand times that of 
DGEAR. 
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Figure (2.6-4) DGEAR 
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Table (2.6-2) Comparison Table for DVERK and DGEAR 


| DVERK j DGEAR | 

I I I 

| Verner-Runge-Kutta j Adam's Predictor-Corrector | 

j 5th, 6th order Method j Variable order (GEAR's) Method j 


a - 1.0 E + 4 


00.00.30 seconds (CPU) | | 

6.83 E - 6 MGRE * j j 

449 KBS ** j j 


o = 1 . 0 E + 6 


| 00.03.26 seconds (CPU) | 00.00.05 seconds (CPU) 

| • -1.0 E - 2 MGRE * | -1.1 E - 5 MGRE * 

j 231468 KBS ** j 7492 KBS ** 


o = 1.0 E + 14 


I 

I 

I 

I 


! 00.00.06 seconds (CPU) | 
| 2.0 E - 2 MGRE * | 
j 1346 KBS ** | 

I I 


* MGRE is defined, using equations (2.6-3) through 
(2.6-6), as the maximum relative error for either i^(t) or 

i 2 (t). 

** KBS is defined as kilobyte-sec which is equal to 
storage (in kilobytes) times CPU time (in seconds). 
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The true cost of each routine is directly proportional to 
its use of both memory and execution time. The cost is a 
function of the product of storage used and CPU time. This 
cost can be expressed in terms of kilobyte- seconds (KBS) 
utilized by each routine to obtain its numerical 
approximation of the exact solution. For example, DVERK 
required approximately thirty times the KBS used by DGEAR to 
obtain the numerical solution for the case mentioned above. 
Even though DGEAR used about twice as much storage as DVERK, 
in this case, DVERK required about sixty-five times longer 
than DGEAR to execute because of its time step restriction 
mentioned earlier. Consequently, for this case, DGEAR was 
approximately thirty times less expensive and had a MGRE 
three orders of magnitude smaller than DVERK. 

Finally, a case illustrating extreme stiffness, where o = 
14 -3 

1.0 E and TOL = 1.0 E , was examined. The results for 

this case are given in the last row of Table (2.6-1). 

Notice that DVERK was not able to track the exact solution 

for this large stiffness. This case illustrates that DGEAR 

can handle much stiffer systems that DVERK. Thus, even 

though the order of stiffness was increased by a factor of 

14 

eight to o = 1.0 E , DGEAR still tracked the exact solution 
as well as DVERK did for the less stiff case when 
o = 1.0 E 6 . Furthermore for this case, DGEAR also excelled 
in the other two categories (of CPU seconds of execution 
time and KBS) even when the stiffness was increased by eight 
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orders of magnitude, as shown in Table (2.6-2). 

The comparison of the two IMSL routines, DVERK and DGEAR, 
yields the conclusion that DGEAR is far superior to DVERK 
when solving stiff systems of differential equations. Even 
though DGEAR uses more memory than DVERK, DVERK' s excessive 
CPU time requirements far outweigh this memory advantage. 
Further, DVERK is not only more expensive but is less 
accurate (as summarized in Table (2.6-2)) to use than DGEAR, 
when integrating stiff systems of differential equations. 
Thus, DGEAR was chosen over DVERK as an appropriate 
integration routine ,for solving the differential equations 
associated with electromechanical systems employing solid 
state switching devices. 

DGEAR was used for integrating the state equations for 
all the simulation runs presented later in this 
dissertation. These runs were made economically feasible 
because of the employment of DGEAR as well as the diode 
routine discussed previously within Section 2.1. 
Specifically, the diode algorithm does not interrupt DGEAR' s 
control of the variable integration step size except for 
diode zero crossings. Hence, it allows DGEAR to integrate 
as economically as it can (maximizing its step size) for the 
majority of the simulation duration. This combination is 
most suitable for the types of machine systems considered in 
this report. 
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2.7 Verification of Modeling Approach . 


This section deals with the verification of the modeling 
approach presented in previous sections of this chapter. 
The system chosen for this purpose was an advanced 
electronically-commutated brushless dc motor designed for 
electric-vehicle propulsion, [20, 21, 22]. The performance 
and operating voltage and current waveforms of this system 
are well documented in these references. A block 

diagram of the propulsion system is presented in Figure 
(2.7-1). 

The verification of the modeling approach presented in 
this section is based upon the favorable comparison of 
actual system measurements and simulated results. This 
section presents verification of the modeling approach's 
ability to simulate systems containing transistor and diode 
switching devices as well as mutual inductances. 

A detailed schematic of the electric vehicle propulsion 
unit is given in Figure (2.7-2). The corresponding lumped 
parameter network model is given in Figure (2.7-3). 
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Figure (2.7-1) Block Diagram of the Electric Vehicle Propulsion Unit 
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The chosen test runs correspond to the rated 15 hp 
operation of this system and additional simulation data are 
presented in Table (2.7-1). The diodes and transistors of 
Figure (2.7-2) were all modeled as seperate (nonlinear) 
piece-wise linear resistors as shown in Figure (2.7-3). The 
oscillograms and other test data presented in this section 
were obtained with the brushless dc machine coupled to a 
dynamometer. 

The parameters of the electric-vehicle propulsion unit as 
well as the control signals for all the transistors were 
obtained from reference [21]. The emf voltage waveforms 
were obtained from a field analysis [38] while the 
inductance values (self and mutuals) were obtained from the 
motor. The machine-inductance values (self and mutuals) 
were obtained from test measurements made on the actual 
machine modeled. 

Two operating modes were simulated using the modeling 
approach presented in this report. The first mode simulated 
was motoring without chopping. The second mode was motoring 
with chopping. Details of these two motoring modes can be 
obtained in reference [21], 
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Figure (2.7-2) Detailed Schematic and Corresponding Network 

Model of the Electric Vehicle Propulsion 
Unit 
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Figure (2.7-3) 



>d Parameter Network Model of 
:ric Vehicle Propulsion Unit. 



the 
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Table (2.7-1) Simulation Data for Electric Vehicle 
Propulsion Unit (Verification Runs) 


Brushless DC Motor Components : 

field: Samarium Cobalt permanent magnet effects are 

reflected as open circuit emf in armature. 

armature: Three-phase wye-connected with floating neutral 

phase to neutral winding resistance = 6.8 milli-ohms 

phase to neutral winding inductance values: 

self = 47.82 micro-henries 

mutual = 4.18 micro-henries 

phase to neutral winding open circuit emf functional 

representation is a truncated fourier series with 

coefficients for the sine terms (a ) and the 

n 

cosine terms (b n ), for the n-th harmonic, given 
below: 

n a /RMVEL b /RMVEL 

n n 


1 

74.400 

3.550 

3 

-1.420 

-0.221 

5 

0.000 

0.000 

7 

0.553 

-0.619 

9 

0.664 

-0.278 

11 

0.373 

-0.141 

13 

-0.037 

-0.014 


where: 


RMVEL = (n RpM /60.0) * 2 * pi * 10 


Power Conditioner Components : 


diode: 


"on" resistance = 6.4 milli-ohms 

"off" resistance = 20.0 kilo-ohms 


transistor: 


"on" resistance 
"off" resistance 


6.4 milli-ohms 
20.0 kilo-ohms 


choke: winding resistance 

self inductance 


24.0 milli-ohms 
1.5 milli-henries 


battery: 


thevenin equivalent resistance = 5.0 milli-ohms 

thevenin equivalent voltage = 120.0 volts 
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The first motoring mode was simulated at the rated 15 hp 
operating point. For this mode, the chopper transistor Q 7 
is latched-on, in the saturation mode, andQg is latched in 
its cutoff mode. Figure (2.7-2). This mode results in the 
full battery voltage being applied continuously across the 
inverter/converter bridge. The motoring results from 
periodic application of the battery voltage across the 
machine's terminals (line to line) in a specific sequence 
that is controlled by transistors through Qg. The 
oscillogram of the phase current of the electric-vehicle 
motor given in Figure (2.7-4) matches very closely the 
simulated waveform shown in Figure (2.7-5). The model was 
used also to predict the voltages across the transistor and 
diode switches. The oscillogram. Figure (2.7-6), and the 
predicted collector to emitter voltage waveform. Figure 
(2.7-7), of one of the inverter transistors are given for 
comparison. Inspection of these two waveforms reveals very 
close agreement between measured and simulated results for 
this first motoring mode of operation. 

Another motoring mode was simulated and the results are 
compared with available oscillograms. This mode is called 
chopping. What is meant by chopping is that the conduction 
status of transistors and Qg are controlled by logic base 
control signals which are determined by monitoring the 
magnitude of the current labeled CM, in Figure (2.7-3). 
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BRANCH 1. X BASE IN SEC. =0.0100 


Figure (2.7-5) Predicted Waveform of the Phase Current 
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Figure (2.7-7) Predicted Waveform of the Inverter Switch 

Voltage 



This feedback control is accomplished by external circuitry 
for the actual system and by a logical subroutine within the 
program. Both of the simulated waveforms and their 
corresponding oscillograms presented in Figure (2.7-8) 
through (2.7-11), reveal excellent correlation. Hence, both 
simulated motoring modes (with and without chopping) agree 
excellently with the actual system measurements. 

Both of these motoring modes (with and without chopping) 
as well as additional verification of a portion of this 
modeling approach was presented earlier in reference [20] 
for the same system and will not be repeated here. In all 
cases, the agreement between measured and simulated results 
was excellent. What was not shown in reference [20] was the 
added feature to have a capability for handling mutual 
inductive coupling between the tree and cotree of a network. 
Also, the incorporation of an implicit versus an explicit 
integration routine , [19], reveals excellent correlation in 
every branch, with the improvement of less computation time. 
Also, the employment of the diode algorithm presented in 
Section 2.5 along with the implicit routine account for the 
improvement in computational efficiency and accuracy. The 
verification of the capability to handle mutual inductances 
between the tree' and cotree of the overall system model will 
be discussed next. 


130 





Figure (2.7-9) Predicted Waveform of the Phase Current 
(Motoring with Chopping Control) 
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Figure (2.7-11) Predicted Waveform of the Inverter Switch 

Voltage (Motoring with Chopping Control) _ 


n 
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Further verification of the proposed modeling approach 
shall address the equality of representing the balanced 
machine’s inductances as a lumped diagonal inductive matrix 
versus a full symmetric inductive matrix. Comparison of the 
simulated results obtained by both inductive matrices 
indicates the validity of the proposed matrix algebra. This 
novel approach includes the added feature which can 
automatically account for mutual coupling between the tree 
and cotree of the system network model, for any arbitrary 
windings. Further, the proposed matrix algebra is not 
restricted to balanced machine operation and is therefore 
directly amendable to fault analysis. Figures (2.7-5) and 
(2.7-7) represent the simulated phase current and inverter 
switch voltage, respectively, when the mutual terms are 
lumped with the self inductances. 

Figures (2.7-12) and (2.7-13) represent the simulated phase 
current and inverter switch voltage, respectively, when the 
inductance matrix is full. A comparison of these waveforms 
reveals excellent correlation between these two sets of 
results. Details of the algebraic manipulations relating to 
this discussion can be seen in equations (2.3-62) and 
(2.3-64). The mutual inductance terms (those between the 
tree and cotree) are coupled via the hybrid matrix pre- and 
postmultiplications as shown in equation (2.3-62). 




Figure (2.7-12) Predicted Waveforms of the Phase Current 




Figure (2.7-13) Predicted Waveforms of the Inverter Switch 
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To assess the validity of this verification, the 
influence of the mutuals upon the overall system performance 
was determined by simulating this same system while 
completely eliminating the mutuals from the full inductance 
matrix. The simulated waveforms, given in Figure (2.7-14), 
indicate an approximate 3 percent increase in the magnitude 
of the phase current for the case when the mutual terms were 
removed. The effects of the mutuals are mild, but 
noticeable, in this case and should not in general be 
ignored when simulating such systems considered in this 
report. The physical reason for the relatively weak (3 
percent) effect of the mutuals is due to the relatively 
large effective air gap (which magnetically couples each of 
the stator windings). The effective air gap is larger 
because the permeability of the samarium-cobalt permanent 
magnets, located on the rotor, is very close to that of air. 
Hence, the mutual coupling between any of the armature 
windings is weak as indicated by these simulated results. 

The next point of interest is how much influence the 
"off" resistance values of the diodes have upon the 
simulated waveshapes and magnitudes as well as the effect 
upon power delivered from the battery power source. Three 
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BRANCH 1 . X BASE IN SEC. =0.0100 


Figure (2.7-14) Predicted Waveforms of the Phase Current 

(Mutuals Removed) 
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cases were run for different "off" resistance values but 
only the two extreme (minimum and maximum) cases are 
displayed. The low value of "off" resistance was set equal 
to 200 ohms and the maximum value was set at 20,000 ohms. 
This hundred fold difference in "off" resistance values 
resulted in less diode leakage current during the "off” 
state but the waveshapes were effectively unaltered. These 
runs were made by only varying the "off” resistance of the 
diodes not the transistors. 


The two sets of diode waveforms agree closely for these 
two cases, as shown in Figures (2.7-15) through (2.7-18). 
The first set of diodes waveforms are for branch 23, the 
positive-Bus phase-A diode. The second set of voltage and 
current waveforms are for branch 22, the flyback diode. The 
actual location of these two branches in the overall network 
can be seen in Figures (2.7-2) and (2.7-3). Both sets of 
waveforms show little alteration in magnitude and waveshape 
as a result of the two values of This is also true 
for the next set of figures to be discussed. The average 
battery current can be obtained from Figures (2.7-19) and 
(2.7-20) for these two values of R off .. The ripple peak to 
peak is unaltered, hence the difference in the YBASE values 
of these Figures’ is equal to the difference in their average 


currents. 




BRANCH 23. X BASE IN SEC. =0.0100 


Figure (2.7-15) Predicted phase A (+ bus) diode waveforms 

for R Q £f = 20.0 kilo-ohms 
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Figure (2.7-17) 


22 . 


X BASE IN SEC. =0.0100 


Predicted (flyback) diode waveforms for 
R = 20.0 kilo-ohms 




BRANCH 22 . X BASE IN SEC. =0.0100 


Figure (2.7-18) 


Predicted (flyback) diode waveforms for 

R = 0.20 kilo-ohms 
off 
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BRANCH 4. X BASE IN SEC. =0.0100 


Figure (2.7-19) Predicted battery waveforms for R __ = 20.0 

off 

kilo-ohms 
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In conclusion, the modeling approach presented here gives 

excellent results in comparison with test data for the types 

of systems considered in this report. It was shown that the 

voltage and current waveshapes are not affected noticeably 

by the choice of the ,, off ,, resistance, above 200 ohms in 

this case, for the solid state switches. However it was 

demonstrated that the ''off” resistance is significant when 

calculating system efficiencies. Therefore, if efficiency 

calculations are not of interest when simulating such 
% 

systems, the value of can be relaxed. This relaxation 

is desirable since the computation cost increases with the 
increase of the value chosen for R Qff . Furthermore, the 
ability to include mutual inductances in the machine model 
was demonstrated and verified against test data. 
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3.0 APPLICATION OF THE ACTUATOR MODEL TO 
VARIOUS MACHINE-POWER CONDITIONER CONFIGURATIONS 


The actuator model developed in the previous chapter will 
now be applied to four additional PSA-machine 
configurations. These are a delta, open-delta, a wye and a 
wye with a shorted turn. It is assumed here that each of 
these systems has a separate chopper network. In general, 
however, the transistor switching logic routine can handle a 
number of different system configurations. The various 
options available are listed in Table (3.0-1). These 
options increase the flexibility of this model by allowing a 
large variety of possible configurations with a minimum 
amount of USER supplied input data. 

In order to properly use this model, several conventions 
must be followed, see Figure (3.0-1). First the 
relationship between the phase and line currents for the 
delta type machines must satisfy the following relationship 
with respect to the line currents i^, anc * igJ 



these constraints result from 
that machine. 


(3.0-1) 

(3.0-2) 

(3.0-3) 

KCL taken at the terminals of 
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Table (3.0-1) PSA Options 


N. ' Variable 
N. Value 

Variable's. 
Name . 


1 

2 

3 





■ 

INVTP 


Wye 

Connected 

Machine 

Delta 

Connected 

Machine 

Open-Delta 

Connected 

Machine 

I CHOP . 


Separate 

Chopper 

Chopping 
Performed 
by Inverter 

- 

ICHOPT 

1 

Hysteresis 

PWM 

- 

MODE 

1 

Motoring 

Re gen 

Plugging 
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ll ia 



WYE 

CONNECTION 



ll ia 



*2 ib 

— • > • 


OPEN DELTA 
CONNECTION 



Figure (3.0-1) Assumed Relationships Between Phase and Line 

Currents 
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For the wye connected machine, the phase currents are 
equal to the line currents. This situation is shown in 
Figure (3.0-1). It is also assumed that the fundamental 
components of the forcing functions obey the following phase 
relationships: 


e a = E m sin < 5t + V 


(3.0-4) 

e b = E m sin < 5t + 9 X " 

120) 

(3.0-5) 

e = E sin( St + 0 . - 
cm x 

240) 

(3.0-6) 

The phase voltages 

must be of this 

form otherwise the 


transistor switching signals would be incorrect. 

Another factor which must be taken into account when 
preparing a simulation run is the assignment of labels for 
the power transistors, Q^, Q 2 ,...,Q^ 2 . This information is 
given in Table (3.0-2). Notice that each of the inverter 
transistors is associated with one of the line currents i^, 
i 2 , or ig. The sign of the line current into the machine, 

( + ) or (-), is determined by which of the two transistors, 
attached at each machine terminal, is "on". 

Similarily, the assumed numbering of the chopper transistors is 
given in Table (3*0-3) • Transistors Q7 or Q13 are used for motoring 
while Q8 or Q14* are used for regeneration. 

The line currents, in all cases, are assumed to be rectangular blocks 
of 120 electrical degrees duration. The line current i^ is initiated 30 
degrees after the emf of phase (A) passes through zero in the positive 
direction. The second and third line currents are initiated in a 
similar fashion 120 and 240 degrees later, respectively. 
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TABLE (3.0-2) Inverter Trans is tor Numbering 


MACHINE/ 

INVERTER 

TYPE 


ACTIVE TRANS 1ST 

ORS 

SIGN OF INJECTED 
LINE CURRENT 

mm 

1 2 

2 3 



Ql 

Q2 

Q3 

( + ) 

WYE 


* 






Q4 

Q5 

Q6 

(-) 

DELTA 


Ql 

Q2 

Q3 

( + ) 



Q4 

Q5 

Q6 

(-) 



Ql 

Q3 

Q5 

( + ) 



Q2 

Q4 

Q6 

(-) 

OPEN DELTA 








Q7 

Q9 

Qll 

(-) 



Q8 

Qio 

Q12 

( + ) 
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TABLE (3*0-3) Chopper Transistor Numbering 


MB 

CHOPPER TRANSISTORS 

WITH SEPARATE CHOPPER 

WITHOUT SEPARATE CHOPPER 

WYE 

07,08 

Q1 -06 

DELTA 

07,08 

Q1-Q6 

OPEN 

DELTA 

Q1 3,014 

Q1-Q12 
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The parameters used in the four simulations which are presented 


next are based upon the system analyzed in the previous chapter. 

In the case of the delta connected machine, the number of turns 
per phase is increased by a factor of •TT. It was assumed that the 
copper volume remained constant, therefore one can write 

(3.0-7) 
(3.0-8) 


Vol CU = *y A y ~ *4*A 
^‘y = ‘a 


where l , A y , £^, are the lengths and crossectional areas 
of wye and delta machines, respectively. Using this 
information one can relate the winding resistances of the 
two cases as follows: 


R a = 3R y (3.0-9) 

The inductances also increase by the square of the turns 
ratio, that is 

L a = 3L y (3.0-10) 

Once all of the machine parameters were scaled properly, 
these parameters were then substituted into the model and 
the simulated waveforms were obtained. The network graphs 
corresponding to these three cases are given in Figures 
(3.0-2), (3.0-3) and (3.0-4). Representative samples of 

simulated branch voltage and current waveforms are given in 
Sections 3.1 through 3.3 for the network graphs (models) of 
Figures (3.0-2) through (3.0-4), respectively. Saturation 
effects in the machine were neglected in all three cases. It was also 
assumed that there were no faults . The machine speed in all three 
cases was 7750 rpm (4 poles) . 

The effects of a shorted turn in phase (A) of the wye connected 
machine are given in Section 3*4. 
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Figure (3.0-2) 


r » 12 . 



Network Graph of WYE Configuration 
(Unfaulted) 
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Figure (3.0-4) Network Graph of Open-DELTA Configuration 

(Unfaulted) 








3.1 Simulated Waveforms of the WYE Configuration 

The simplified (with chopper) schematic diagram of 
Figure (1.0-1) for the WYE configuration, and its 
corresponding complete network graph (including chopper) of 
Figure (3.0-2) should be referred to in reviewing the 
results of simulations presented in this section. These 
simulations represent unfaulted normal motoring operation. 
The input data used in obtaining these results are given in 
Table (3.1-1). Identification of all the displayed voltage 
and current simulated waveforms is given in Table (3.1-2). 

The calculated average electromagnetic machine power during the 
last cycle is 12, 6l? watts. The corresponding copper losses in the 
machine totaled 218.7 watts. 
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TABLE (3.1-1) Input Data for Network Model of the WYE 


Configuration (Unfaulted) 



Three Phase Balanced Machine Model 


Branch #12 Phase(A) Winding Resistance=6. 8 milli-ohms 
Branch #11 Phase (A) Winding Inductance=52 micro-henries 


Inverter/Converter Bridge & 


Switching Components: 


Branch #10 (+Bus) Phase (A) Transistor 
Saturation Resistance "on" = 6.4 milli-ohms 
Cutoff Resistance "off" =2.0 killo-ohms 

Branch #23 (+Bus) Phase (A) Diode 
Forward-Biase Resistance "on" = 6.4 milli-ohms 
Reverse-Biase Resistance "off"= 2.0 killo-ohms 

Choke : 


Branch #15 Choke Resistance 
Branch #17 Choke Inductance 

DC Power Source: 


24 milli-ohms 
1.5 milli-henries 


Branch #5 Filter Capacitor = 7.8 milli-f arrads 
Branch #4 Battery ( Internal) Voltage = 120.0 Volts 
Branch #16 Battery ( Internal ) Resistance = 5 milli-ohms 




Table (3*1-2) Branch Identification WYE Configuration 
( Unfaulted) 


Figure 

Branch 

Branch Identification 

Number 

Number 



HH 



3*1-1 

1 

Phase (A) Back Emf 

3*1-2 

11 

Phase (A) Winding Inductance 

3*1-3 

10 

Phase (A) Transistor (+Bus) 

3*1-4 

23 

Phase (A) Diode (+Bus) 

3*1-5 

7 

Phase (A) Transistor (-Bus) 

3*1-6 

29 

Phase (A) Diode (-Bus) 

3*1-7 

4 

Battery Internal Voltage 

3*1-8 


Chopper Inductance 

3*1-9 


Chopper (+Bus, 

Transistor 

3*1-10 


Chopper (+Bus) 

Diode 

3*1-11 

«Dn 

Chopper (-Bus] 

Transistor 

3*1-12 

28 

Chopper (-Bus) 

Diode 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AHPS 




Figure (3.1-1) Simulated Current and Voltage Waveforms for 

Branch #1, WYE Configuration (Unfaulted) 
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BRANCH VOLTAGE IN VOLTS . BRANCH CURRENT IN RHPS 
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, YBRSE- 1.47 xlQ 3 XBASE- 4.00x10"* 




Figure (3.1-3) Simulated Current and Voltage Waveforms for 

Branch #10, WYE Configuration (Unfaulted) 
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line IN SECONDS 


Figure (3.1-4) Simulated Current and Voltage Waveforms for 

Branch #23, WYE Configuration (Unfaulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AHPS 




Figure (3.1-5) Simulated Current and Voltage Waveforms for 

Branch #7, WYE Configuration (Unfaulted) 
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in mps 







. YBRSE- 1. 49x10’ XBASE- 4.00x10"’ 

1 I 




Figure (3.1-7) Simulated Current and Voltage Waveforms for 

Branch #4, WYE Configuration (Unfaulted) 
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BRANCH VOLTAQE IN V0LT3 , BRANCH CURRENT IN RHP3 




Figure (3.1-8) Simulated Current and Voltage Waveforms for 

Branch #17, WYE Configuration (Unfaulted) 
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Figure (3.1-9) Simulated Current and Voltage Waveforms for 

Branch #20, WYE Configuration (Unfaulted) 
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BRUNCH VOLTRQE IN VOLTS , BRANCH CURRENT IN RHPS 




Figure (3.1-10) Simulated Current and Voltage Waveforms for 

Branch #21, WYE Configuration (Unfaulted) 
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BRANCH VOLTHQE IN VOLTS . BRANCH CURRENT IN AMPS 


YBflSE- S.SSxlQ XBASE- 4.00x10' 



[ij 


Tire IN SECONDS 


YBflSE- 1. 19x10* XBASE- 4.00x10"* 



Tire INSBOOCS 


Figure (3.1-11) Simulated Current and Voltage Waveforms for 

Branch #6, WYE Configuration (Unfaulted) 
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BRANCH VOLTAQE IN VOLTS , BRANCH CURRENT IN AMPS 




Figure (3.1-12) Simulated Current and Voltage Waveforms for 

Branch #28, WYE Configuration (Unfaulted) 
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3.2 Simulated Waveforms of the DELTA Configuration 

The simplified (with chopper) schematic diagram of 

Figure (1.0-2) for the DELTA configuration, and its 
corresponding complete network graph (including chopper) of 
Figure (3.0-3) should be referred to in reviewing the 
results of simulations presented in this section. These 
simulations represent unfaulted normal motoring operation. 
The input data used in obtaining these results are given in 
Table (3.2-1). Identification of all the displayed voltage 
and current simulated waveforms is given in Table (3.1-2). 

The calculated average electromagnetic machine power during . the 
last cycle is 12,594 watts. The corresponding copper losses in the 
machine totaled 217.4 watts. These values are almost identical to the 
ones obtained from the WYE connected machine. This is to be expected 
since the machine parameters were chosen such that their terminal per- 
formances are identical . 
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Table (3.2-1) Input Data for Network Model of the DELTA 


Configuration (Unfaulted) 


i 1 

I ’ I 

I Three Phase Balanced Machine Model: j 


Branch #11 Phase(A) Winding Resistance=20. 4 milli-ohms 
Branch #17 Phase (A) Winding Inductance=156 micro-henries 

Inverter/Converter Bridge & Chopper Switching Components : 

Branch #7 (+Bus) Phase(A) Transistor 
Saturation Resistance "on" = 6.4 milli-ohms 
Cutoff Resistance ’'off*' =2.0 killo-ohms 

Branch #22 (+Bus) Phase(A) Diode 
Forward-Biase Resistance "on” = 6.4 milli-ohms 
Reverse-Biase Resistance "off"= 2.0 killo-ohms 

« 

Chopper Choke : 


Branch #14 Choke Resistance = 24 milli-ohms 
Branch #16 Choke Inductance = 1.5 milli-henries 

DC Power Source : 

Branch #5 Filter Capacitor = 7.8 mi Hi- farads 
Branch #4 Battery ( Internal) Voltage = 120.0 Volts 
Branch #15 Battery( Internal ) Resistance = 5 milli-ohms 


i i 







Table (3.2-2) Branch Identification DELTA Configuration 

(Unfaulted) 


Figure 

Number 


Branch 

Number 


Branch Identification 


3.2-1 | 

1 

| Phase (Aj 

Back Emf 

3.2-2 j 

17 

| Phase (a' 

Winding Inductance 

3.2-3 j 

7 

| Phase (A 

Transistor (+Bus) 

3.2-4 | 

22 

| Phase (A ( 

Diode (+Bus) 

3.2-5 j 

8 

| Phase (a' 

Transistor (-Bus) 

3.2-6 j 

28 

| Phase (A; 

Diode (-Bus) 

3.2-7 j 

4 

j Battery 

Internal Voltage 

3.2-8 j 

16 

j Chopper 

Inductance 

3.2-9 j 

31 

1 Chopper 

(+Bus) Transistor 

3.2-10 j 

20 

S Chopper 

(+Bus) Diode 

3.2-11 | 

6 

j Chopper 

(-Bus) Transistor 

3.2-12 j 

27 

j Chopper 

(-Bus) Diode 
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YBfiSE- 1.07*10’ 






BRflNCH V0LTR9E IN VOLTS , BRANCH CURRENT IN flHPS 




Figure (3.2-2) Simulated Current and Voltage Waveforms for 

Branch #17, DELTA Configuration (Unfaulted) 


9 




BRANCH VOLTflGE IN VOLTS 


-1 

0 


TIME IN SECONDS 


I 



Figure (3.2-3) Simulated Current and Voltage Waveforms for 

Branch #7, DELTA Configuration (Unfaulted) 
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YQLTR8E IN VOLTS , BRflNCH CURRENT IN flNPS 



TItC IN SECONDS 


Figure (3.2-4) Simulated Current and Voltage Waveforms for 

Branch #22, DELTA Configuration (Unfaulted 
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BRANCH VOLTA* IN V0LT3 , BRANCH CURRENT IN AMPS 




XBASE' 




-I 


TIME IN SECONDS 


Figure (3.2-5) Simulated Current and Voltage Waveforms for 

Branch #8, DELTA Configuration (Unfaulted) 
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BRANCH YOLTRQE IN VOLTS 



Figure (3.2-6) Simulated Current and Voltage Waveforms for 

Branch #28, DELTA Configuration (Unfaulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AMPS 
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. YBflSE- 1.20x10* XBASE- 4.00 xlO* 

I r ~ — — — — 1 
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-l 
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Tire IN SECONDS 


1 


Figure (3.2-7) Simulated Current and Voltage Waveforms for 

Branch #4, DELTA Configuration (Unfaulted) 


184 - 






Figure (3.2-8) Simulated Current and Voltage Waveforms for 

Branch #16, DELTA Configuration (Unfaulted) 
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Figure (3.2-9) Simulated Current and Voltage Waveforms for 

Branch #31, DELTA Configuration (Unfaulted) 
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BRANCH YOLTME IN VOLTS BRANCH CURRENT IN AMPS 


YBflSE- 6.02x10“* XBASE- 4.00x10**' 



Tire IN SGCGNDS 


YBRSE- 1.20xl0 a XBASE- 4.00x10 



TOC IN SECONDS 


Figure (3.2-10) Simulated Current and Voltage Waveforms for 

Branch #20, DELTA Configuration (Unfaulted) 
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Figure (3.2-11) Simulated Current and Voltage Waveforms for 

Branch #6, DELTA Configuration (Unfaulted) 
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BRANCH VOLTA06 IN VOLTS i BRANCH CURRENT IN AMPS 




Figure (3.2-12) Simulated Current and Voltage Waveforms for 

Branch #27, DELTA Configuration (Unfaulted) 
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3.3 Simulated Waveforms of the OPEN-DELTA Configuration 

The simplified (with chopper) schematic diagram of 
Figure (1.0-3) for the DELTA configuration, and its 
corresponding complete network graph (including chopper) of 
Figure (3.0-4) should be referred to in reviewing the 
results of simulations presented in this section. These 
simulations represent unfaulted normal motoring operation. 
The input data used in obtaining these results are given in 
Table (3.3-1). Identification of all the displayed voltage 
and current simulated waveforms is given in Table (3.3-2). 

The calculated average electromagnetic machine power during the 
last cycle is 18,40? watts. The corresponding copper losses In the 
machine totaled 515 watts. 
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Table (3.3-1) Input Data for Network Model of the 


OPEN-DELTA Configuration (Unfaulted) 



Branch #15 Phase (A) Winding Resistance=20.4 milli-ohms 
Branch #20 Phase(A) Winding Inductance=156 micro-henries 



Branch #7 (+Bus) Phase(A) Transistor 
Saturation Resistance "on" = 6.4 milli-ohms 
Cutoff Resistance "off" = 2.0 killo-ohms 


. Branch #26 (+Bus) Phase (A) Diode 

Forward-Biase Resistance "on" =6.4 milli-ohms 
Reverse-Biase Resistance. "off"= 2.0 killo-ohms 


Chopper Choke : 

Branch #14 Choke Resistance = 24 milli-ohms 
Branch #19 Choke Inductance = 1.5 milli-henries 

DC Power Source : 


Branch #5 Filter Capacitor = 7.8 milli-farads 1 
Branch #4 Battery ( Internal) Voltage = 120.0 Volts | 
Branch #18 Battery ( Internal ) Resistance = 5 milli-ohms | 
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Table (3.3-2) Branch Identification OPEN-DELTA Configuration 

(Unfaulted) 


Figure 

Number 

Branch 

Number 

Branch Identification 




3.3-1 

1 

Phase (A) Armature EMF and Phase Current 

3.3-2 

20 

Phase (A) Inductive Voltage and Phase Current 

3.3-3 

7 

Phase (A) (+ Bus) Transistor 

3.3-4 

26 

Phase (A) (+ Bus) Diode 

3.3-5 

9 

Phase (A) (-Bus) Transistor 

3.3-6 

39 

Phase (A) (-Bus) Diode 

3.3-7 

4 

Battery Internal Voltage and Current 

3.3-8 

19 

Chopper Inductive Voltage and Current 

3.3-9 

23 

Chopper (+ Bus) Transistor 

3.3-10 

24 

Chopper (+ Bus) Diode 

3.3-11 

6 

Chopper (- Bus) Transistor 

3.3-12 

37 

Chopper (- Bus) Diode 








BRANCH V0LTA8E IN VOLTS , BRANCH CURRENT IN ANP3 



Figure (3.3-1) Simulated Current and Voltage Waveforms for 

Branch #1, OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH YOLTftGE IN VOLTS 



Figure (3.3-2) Simulated Current and Voltage Waveforms for 

Branch #20, OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAGE IN VOLTS . BRANCH CURRENT IN ANPS 



YBflSE- 1. 22x10* XBASE- 4.00x10"* 



Figure (3.3-3) Simulated Current and Voltage Waveforms for 

Branch # 7 , OPEN-DELTA Configuration 
(Unfaulted) 
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IN VOLTS 
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TIME IN 9CC0C3 
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. YBflSE- 1.22*10* XBASE- 4.00*10* 

l . — - - — — — 



Eigure (3.3-4) Simulated Current and Voltage Waveforms for 

Branch #26, OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AMPS 




Figure (3.3-5) Simulated Current and Voltage Waveforms for 

Branch #9 , OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH 




Figure (3.3-6) Simulated Current and Voltage Waveforms for 

Branch #39, OPEN-DELTA Configuration 
(Unfaulted) 
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9 



-1 


YBASE- 1.20*10* XBASE* 4.00*10"* 


TIfE W3BDQW3 


Figure (3.3-7) Simulated Current and Voltage Waveforms for 

Branch #4, OPEN-DELTA Configuration 
(Unfaulted) 


■s 
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Figure (3.3-8) Simulated Current and Voltage Waveforms for 

Branch # 19 , OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAGE IN VOLT3 , BRANCH CURRENT IN AMPS 



. XBASE- 1.20x10* XBASE- 4.00 xlO* 

1 I 


0 


TIME IN SECONDS 


Figure (3.3-9) Simulated Current and Voltage Waveforms for 

Branch #23/ OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAQE IN VOLTS , BRANCH CURRENT IN AMPS 
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. YBRSE- 1.20*10* XBASE- 4.00x10^ 
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Tire IN SECONDS 


1 


Figure (3.3-10) Simulated Current and Voltage Waveforms for 

Branch #2.4, OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VQLTRQE IN VOLTS , BRANCH CURRENT IN RHP8 


YBflSE" 5.97x10* XBRSE- 4.00*10°* 
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. YBflSE- 1. 19x10’ XBASE- 4.00x10'* 

1 r- — — 1 
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-1 


0 


Tire in scores 
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Figure (3.3-11) Simulated Current and Voltage Waveforms for 

Branch #6, OPEN-DELTA Configuration 
(Unfaulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AHP9 
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Figure (3.3-12) Simulated Current and Voltage Waveforms for 

Branch #37, OPEN-DELTA Configuration 
(Unfaulted) 
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3.4 Simulated Waveforms of the WYE Configuration 
( Containing Shorted Turn ) 

A simulation of the effects of a shorted turn in a WYE 
connected machine is presented in this section. This 
simulation is otherwise identical to the case presented in 
Section (3.3-1). The corresponding network graph is given 
in Figure (3.4-1) and should be referred to in reviewing the 
results of simulations presented in this section. The 
inductances used to simulate the effects of the shorted turn 
are given in Table (3.4-3). 

The input data used in obtaining these results are given 
« . 

in Table (3.4-1). Identification of all the displayed 
voltage and current simulated waveforms is given in Table 
(3.4-2) . 

The calculated average electromagnetic machine power during the last 
cycle is 12,520 watts. This is only slightly less than that for the un- 
faulted case. The corresponding copper losses, however, are much greater 
(597.15 watts). The losses in the shorted turn are 40 times higher 
(366 watts) than the unfaulted losses (9.1 watts). It must be emphasized 
that saturation effects and changes in winding resistance due to increased 
temperatures were neglected in this analysis. These factors should be 
included in any future efforts in this area. 
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Table (3.4-1) Input Data for Network Model of the WYE 
Configuration ( Faulted ) 


r — i 

I I 

I Three Phase UnBalanced Machine Model: I 


Branch #14 Healthy Portion Phase (A) 7 TURNS 

Winding Resistance* 5,95 milli-ohms 
Branch #13 Shorted Portion of Phase (A) 1 TURN 
Winding Resistance= O .85 milli-ohms 

Branch #15 Phase (B) Winding Resistance= 6.8 milli-ohms 
Branch #16 Phase (C) Winding Resistance= 6.8 milli-ohms 

See Inductance Matrix (Unbalanced Case) Attached in 
Figure (3.4-2) 

Inverter/Converter Bridge & Chopper Switching Components : 

Branch #11 (+Bus) Phase (A) Transistor 

Saturation Resistance "on” = 6.4 milli-ohms 

Cutoff Resistance "off” =2.0 killo-ohms 

Branch #26 (+Bus) Phase (A) Diode 

Forward-Bias Resistance "on” = 6.4 milli-ohms 
Reverse-Bias Resistance "off"= 2.0 killo-ohms 

Chopper Choke: 

Branch #17 Choke Resistance = 24 milli-ohms 
Branch #19 Choke Inductance = 1.5 milli-henries 

DC Power Source : 

Branch #6 Filter Capacitor = 7.8 milli-farads 
Branch #5 Battery ( Internal ) Voltage = 120.0 Volts 
Branch #18 Battery ( Internal ) Resistance = 5 milli-ohms 
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Table (3.4-2) Branch Identification WYE Configuration 

(Containing Short Circuit) 


Figure 

Number 

Branch 

Number 

Branch Identification 

■■I 



3.4-2 

2 

Phase (B) Armature EMF and Phase Current 

3.4-3 

21 

Phase (B) Inductive Voltage and Phase Current 

3.4-4 

27 

Phase (B) (+ Bus) Transistor 

3.4-5 

28 

Phase (B) (+ Bus) Diode 

3.4-6 

9 

Phase (B) (- Bus) Transistor 

3.4-7 

33 

Phase (B) (r Bus) Diode 

3.4-8 

5 

Battery Internal Voltage and Current 

3.4-9 

19 

Chopper Inductive Voltage and Current 

3.4-10 

23 

Chopper (+ Bus) Transistor 

3.4-11 

24 

Chopper (+ Bus) Diode 

3.4-12 

7 

Chopper (- Bus) Transistor 

3.4-13 

31 

Chopper (- Bus) Transistor 

3.4-14 

1 

(H) Phase (A) Armature EMF and Phase Current 

3.4-15 


(H) Phase (A) Inductive Voltage and Phase Current 

3.4-16 


(S) Phase (A) Armature EMF and Phase Current 

3.4-17 


(S) Phase (A) Inductive Voltage and Phase Current 
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BRANCH VOLTAGE IN VOLTS BRANCH CURRENT IN AMPS 



YBRSE- 6. 15x10 1 XBASE- 4.00x10'* 



Eigure (3.4-2) Simulated Current and Voltage Waveforms for 

Banch #2, WYE Configuration (faulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AMPS 



0 
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Figure (3.4-3) Simulated Current and Voltage Waveforms for 

Banch #21, WYE Configuration (faulted) 
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Figure (3.4-4) Simulated Current and Voltage Waveforms for 

Branch #27, WYE Configuration (faulted) 
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Figure (3.4-5) Simulated Current and Voltage Waveforms for 

Branch #28, WYE Configuration (faulted) 
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Figure (3.4-6) Simulated Current and Voltage Waveforms for 

Branch #9, WYE Configuration (faulted) 
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BRANCH VQLTAQE IN VOLTS . BRANCH CURRENT IN RHPS 




Figure (3.4-7) Simulated Current and Voltage Waveforms for 

Branch #33, WYE Configuration (faulted) 
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YBflSE- 1.20x10 
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Tire IN SECONDS 


Figure (3.4-8) Simulated Current and Voltage Waveforms for 

Branch #5, WYE Configuration (faulted) 



BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AHPS 
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YBflSE- 1.26x10* XBASE- 4.00x10 



Eigure (3.4-9) Simulated Current and Voltage Waveforms for 

Branch #19, WYE Configuration (faulted) 
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Figure (3.4-10) 


Branch 


XBASE- 4.00 Kid* 9 


life IN SOCKS 

XBASE- 4.00 xlO* 



TOC IN SZXNS 


Current and Voltage Waveforms for 
, WYE Configuration (faulted) 
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Figure (3.4-11) Simulated Current and Voltage Waveforms for 

Branch #24, WYE Configuration (faulted) 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN ANPS 


. YBflSE- 1.47*10 J XBRSE- 4.00 xlO"’ 
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Figure (3.4-13) Simulated Current and Voltage Waveforms for 

Branch #31, WYE Configuration (faulted) 
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YBflSE- 1.52xlQ 3 XBASE- 4.00x10* 
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BRANCH VOLTAGE IN VOLTS , BRANCH CURRENT IN AHP9 




Figure (3.4-15) Simulated Current and Voltage Waveforms for 

Branch #12, WYE Configuration (faulted) 
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Figure (3.4-17) Simulated Current and Voltage Waveforms for 

Branch #20, WYE Configuration (faulted) 
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4.0 CONCLUSIONS AND RECOMMENDATIONS 


A network topology based modeling approach designed for 
the simulation of electromechanical drive and actuator 
systems has been presented and verified against actual test 
measurements. The salient features of this model are: 

1. The state equations are automatically generated from 
the USER specified network topology using standard 
network graph theory concepts. 

2. The type of control unit used to operate the power 
conditioner can be USER specified from a number of 
preprogrammed controllers. Additional control 
schemes can be implemented, if required, by modifying 
SUBROUTINE QLOGIC. 

3. The program can handle aof number fault scenarios 
including shorted turns and failures in the power 
electronics . 

4. The program can handle mutual inductances and can be 
adopted _to include nonlinearities in the machine 
parameters. 

5. The integration routine employed by this program is 

1 

the IMSL SUBROUTINE DGEAR which is widely available 
and well suited to solving the ill conditioned 
systems of differential equations resulting from 
these types of machine systems. 
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The model was shown to be successful in simulating the 
machines with connections such as WYE, DELTA, OPEN-DELTA, 
and WYE with SHORTED TURN IN ARMATURE WINDING. The model 
was used successfully in determining currents and voltages 
in and across the various branches, respectively. The 
present computer algorithm can be easily supplied with a 
post processor (external subroutines) which would calculate 
input power, output power, and losses in the systems studied 
here. However, the present resources did not permit the 
development of this post processor, which must be added in 
future .research efforts on such machine system simulation 
models. 

There are a number of improvements which would greatly 
expand the capabilities of this model. These include: 

1. A more detailed switch model may be necessary in 

future investigations especially if power conditioner 
efficiencies are of great importance. Such 
improvements may include a more accurate 

representation of the static I-V characteristics as 
well as accounting for switching and base drive 
losses . 

2. A model of SCR switches should also be implemented to 
increase the versatility of the model. 

3. The program should be modified to accept dependent 
sources . 
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4. This model should be linked to a finite element 
machine model when studying the effects of shorted 
turns. This is necessary due to the heavy saturation 
produced by the fault currents. 

5. The program should be streamlined and upgraded to 
enable the USER to perform simulations with a minimum 
amount of input data. This should include a 
preprocessor to automatically generate a proper 
network graph for the system under study. 
Postprocessors to calculate losses, efficiencies, 
^tc. and well as plotting routines to display the 
results would also be useful. 

The results and experimental verification presented in this 
report demonstrate the usefulness and accuracy of this model 
in predicting both the faulted and unfaulted performance of 
electromechanical actuators. Much work remains to be done 
to link this network model with the finite element based 
machine parameter estimation model required for accurate 
simulation of the shorted turn behavior of such systems. 
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